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Effect  of  Particle  Shape  on  the  Mechanical  Behavior  of  Granular  Materials: 

A  Discrete  Element  Study 


John  Ting,  Pricipal  Investigator 
University  of  Massachusetts  Lowell 


EXECUTIVE  SUMMARY 

This  Report  presents  the  results  of  an  Air  Force  Office  of  Scientific  Research-funded  study  on  the  effect 
of  particle  shape  on  the  overall  material  behavior  of  granular  materials  through  a  Discrete  Element  numerical 
study  using  ellipse-shaped  particles.  The  Discrete  Element  Method  (DEM)  possesses  unique  capabilities 
for  studying  the  behavior  of  particulate  materials.  Because  it  is  a  numerical  simulation,  internal  contact 
forces,  particle  orientations  and  average  stresses  and  strains  may  be  monitored  inobstrusively.  Because  the 
DEM  models  each  particle  explicitly,  contacts  may  be  broken  and  renewed,  allowing  for  large  local 
deformations.  However,  past  attempts  at  using  DEM  with  circular  particles  have  resulted  in  simulated 
mechanical  behavior  which  is  considerably  weaker  than  expected  for  natural  granular  materials.  The  first 
portion  of  this  study  involved  calibrating  a  recently  developed  two-dimensional  ellipse-based  computational 
model  (Ting  and  Corkum,  1993)  with  physical  test  data  by  Konishi,  Oda  and  Nemat-Nasser  (1983)  on  oval 
cross-section  rods  and  by  Chapuis  (1976)  on  circular  rods. 

The  ellipse-based  DEM  model  yields  behavior  results  that  are  quantitatively  closer  to  real  soil  than  those 
yielded  by  DEM  models  using  round  particles.  Average  angles  of  internal  friction  ranged  from  26  to  greater 
than  50 degrees,  depending  on  particle  aspect  ratio  and  nature  of  the  anisotropy.  Good  quantitative 
agreement  was  found  with  the  circular  rod  biaxial  shear  data  of  Chapuis  (1976).  Good  qualitative,  but  not 
quantitative,  agreement  was  found  with  the  oval  rod  biaxial  shear  data  of  Konishi  et  al  (1983),  despite 
considerable  effort  in  varying  contact  paramenters  in  the  DEM  model.  The  ellipse-based  DEM  yielded 
consistently  softer  and  weaker  response  compared  to  the  physical  oval  rod  data.  These  differences  are 
attributed  to  the  subtle  but  significant  differences  in  the  shape  of  the  particles  used  in  each  case.  It  is 
recommended  that  validation  of  the  ellipse-based  DEM  should  continue  by  comparing  with  additional 
physical  test  results.  Such  a  study  is  currently  underway  using  the  results  of  interfacial  shear  tests  between 
a  rough  interface  an  assemblage  of  carefully  machined  circular  and  elliptical  rods. 

In  the  second  portion  of  this  study,  additional  numerical  simulations  of  biaxial  shear  tests  were 
conducted  to  determine  the  relative  importance  of  particle  aspect  ratio,  anisotropy,  and  contact  friction. 
Numerical  samples  were  formed  similar  to  the  Konish  et  al  rod  tests,  consisting  of  multiple  sized  particles 
of  varying  particle  aspect  ratio  with  their  major  axes  aligned  along  preferred  bedding  planes.  Individual 
particle  aspect  ratio  (major  axis:minor  axis)  ranged  from  1:1  to  3:1 .  The  samples  were  initially  compressed 
to  isotropic  conditions,  then  sheared  in  biaxial  shear. 

Results  from  these  tests  exhibit  distinctly  different  stress-strain-strength  behavior,  as  well  as  different 
controlling  deformational  mechanisms,  as  a  function  of  both  bedding  plane  angle  and  particle  flatness. 
Interparticle  interlocking,  dilatancy  and  inhibition  of  particle  rolling  were  all  observed  and  significantly 
affect  the  overall  mechanical  behavior.  Overall,  the  samples  with  bedding  orientation  normal  to  the  principal 
stress  direction  exhibited  the  highest  shear  resistance,  as  well  as  the  largest  dilatancy.  Increasing  particle 
flatness  usually  increases  overall  shear  strength,  except  when  the  bedding  plane  coincides  with  the  planes 
of  maximum  obliquity.  Systems  with  randomly  oriented  particles  exhibit  lower  strength  than  systems  with 
bedding  unfavorable  to  shear,  but  exhibit  greater  strength  than  systems  with  bedding  favorable  for  shear. 
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Shear  band  formation  was  observed  in  some  bedding  and  particle  flatness  combinations.  In  systems 
composed  of  flat  particles  with  bedding  parallel  to  one  of  the  planes  of  maximum  obliquity,  the  shear  band 
sometimes  formed  in  the  maximum  obliquity  plane  opposite  to  the  bedding  plane.  When  shear  band 
formation  was  observed,  the  deformation  of  the  assemblage  more  closely  resembles  rigid  body  motion  of 
two  opposing  particulate  wedges  with  the  shear  band  forming  the  slip  plane. 

Based  on  analysis  of  the  contact  statistics  during  shearing,  this  study  confirms  that  particle  rolling  is  a 
dominant  deformational  mechanism  for  systems  with  rounded  particles.  For  systems  with  flatter  particles, 
interparticle  slidng  becomes  more  prominent,  with  the  relative  importance  of  particle  motion  due  to 
translation  and  rotation  becoming  virtually  equal  for  flatter  particles.  For  assemblages  with  low  interparticle 
friction,  particle  sliding  is  the  dominant  deformation  mechanism,  while  particle  rolling  is  the  dominant 
deformation  mechanism  in  assemblages  with  high  interparticle  friction. 

The  DEM  simulations  show  that  assemblages  with  preferred  bedding  possess  more  anisotropic 
distribution  of  contact  normals.  As  particle  angularity  increases,  interlocking  increases  producing  a  more 
anisotropic  fabric.  During  the  course  of  deformation,  the  distribution  of  contact  normal  orientations  tend 
to  rotate  toward  the  direction  of  the  major  principal  stress. 

The  ellipse-based  DEM  model  used  for  the  current  study  has  been  shown  to  be  an  effective  tool  for 
realistically  modeling  the  microstructure  and  overall  mechanical  behavior  of  granular  materials.  Additional 
research  is  planned  or  underway  in  several  areas:  validation  with  systems  composed  of  carefully  machined 
circular  and  elliptical  rods  subjected  to  shear;  modifying  the  code  to  allow  for  fluid-solid  interaction; 
modifying  the  code  to  model  multiphase  materials  such  as  frozen  soil  systems;  validating  the  method  for 
large  scale  problems  such  as  shallow  plate  bearing  capacity;  and  application  of  the  technique  to  specific 
problems  in  granular  flow  and  failure  such  as  pressures  during  silo  filling  and  vehicle-soil  interaction. 

Personnel  who  have  worked  on  this  Research  project,  either  in  a  funded  or  unfunded  capacity  include: 
John  M.  Ting  (Principal  Investigator),  Larry  R.  Meachum,  Jeffrey  D.  Rowell,  Wai-Lim  Chin,  Mahmood 
Khwaja,  John  Sanchez  and  Gregory  Mischel. 


TABLE  OF  CONTENTS 


EXECUTIVE  SUMMARY . . . i 

TABLE  OF  CONTENTS . iii 

1.  INTRODUCTION .  1 

2.  BACKGROUND .  2 

2.1  Effect  of  Particle  Shape  on  Strength .  2 

2.2  Effect  of  Particle  Shape  on  Densification .  2 

2.3  Effect  of  Particle  Shape  on  Pore-pressure  Buildup  .  3 

2.4  Effect  of  Particle  Shape  on  Mechanisms  of  Strength  and  Deformation .  3 

3.  RESEARCH  PROGRAM .  5 

4.  DISCRETE  ELEMENT  METHOD  .  6 

4.1  Derivation  of  Ellipse-based  Discrete  Element  Method  Algorithm  .  7 

4.1.1  Ellipse-ellipse  Intersection  . 8 

4.1.2  Ellipse-Ellipse  Incremental  Contact  Velocities .  11 

4.1.3  Ellipse-line  Intersection  and  contact  velocities .  13 

4.2  Derivation  of  Associated  Properties .  16 

4.2.1  Contact  properties  .  16 

4.2.2  Stress  Tensor .  16 

4.2.3  Strain  Tensor .  17 

4.3  Implementation  of  Ellipse-based  DEM  .  18 

4.4  Code  Validation  .  19 

4.5  Interface  Programs  .  20 

5.  VALIDATION . 22 

5.1  The  Hysteretic  Force-Deflection  Algorithm .  22 

5.2  Nonlinear  Force-Deflection  Relationship .  23 

5.3  Chapuis  (1976)  Test  Data .  27 

5.4  Konishi,  Oda  and  Nemat-Nasser  (1983)  Tests .  32 

5.5  Variation  of  Contact  Properties  .  37 

5.5.1  Variation  of  Linear  Contact  Spring  Stiffness  .  37 

5.5.2  Alternate  Contact  Stiffness  Model  .  38 

5.5.3  Variation  of  Contact  Cohesion  .  38 

5.5.4  Variation  of  Interparticle  Friction  Angle,  ({>„ . 39 

5.5.5  Effect  of  Contact  Damping,  c„ .  39 

5.6  Discussion  of  Validation  Testing  and  Konishi  etc/ Data  .  41 

5.7  Comparison  with  Other  DEM  Results  .  47 

5.7.1  Numerical  Results  from  Rothenburg  and  Bathurst  (1991)  .  47 

5.7.2  Numerical  Results  from  Ng  (1992) .  48 


iii 


6.  SAMPLE  FORMATION  AND  TEST  SIMULATIONS  .  51 

6.1  Sample  Generation  .  51 

6.2  Isotropic  Compression  Simulation .  53 

6.3  Biaxial  Compression  Test  Simulation  .  55 

7.  COMPARISON  WITH  KONISHI,ODA  AND  NEMAT-NASSER  RESULTS  .  59 

7.1  Data  Comparison  for  Samples  with  (J)^  =  26°  .  59 

7.2  Data  Comparison  for  Samples  with  (()„  =  52° .  65 

7.3  Discussion  of  Data  Comparison .  72 

8.  RESULTS  OF  DEM  ANALYSES  .  89 

8.1  Influence  of  Particle  Shape  on  Strength .  89 

8.1.1  Effect  of  Particle  Aspect  Ratio  .  89 

8. 1 .2  Effect  of  Aspect  Ratio  and  Particle  Bedding .  89 

8.1.3  Effect  of  the  Initial  Void  Ratio  on  Strength  .  101 

8. 1 .4  Effect  of  a  Random  Orientation  Fabric . 105 

8.2  Influence  of  Particle  Shape  on  Volumetric  Behavior .  107 

8.2.1  Effect  of  Particle  Aspect  Ratio  .  107 

8.2.2  Effect  of  Particle  Aspect  Ratio  and  Orientation  .  108 

8.2.3  Effect  of  Void  Ratio .  113 

8.2.4  Effect  of  an  Isotropic  Orientation  Fabric .  115 

8.2.5  Deformation  Mechanisms .  116 

8.3  Influence  of  Particle  Shape  on  the  Coordination  Number .  118 

8.4  Effect  of  Sample  Formation .  123 

8.5  Effect  of  Confining  Stress  Oj  .  125 

9.  DISCUSSION  OF  DEM  RESULTS .  129 

9.1  Stress-dilatancy . 129 

9.2  The  Failure  Plane  .  132 

9.3  Contact  Normal  Distribution  .  142 

10.  CONCLUSIONS  AND  RECOMMENDATIONS  .  152 

APPENDIX  I.  REFERENCES  .  155 

APPENDIX  II.  SYMBOLS  .  159 


IV 


1.  INTRODUCTION 


Particle  shape  is  one  of  the  major  factors  determining  the  overall  mechanical  behavior  of  cohesionless 
particulate  matoials.  In  materials  possessing  more  rounded  particles,  the  mechanism  of  rolling  between  particles 
dominates  compared  with  other  deformational  mechanisms  such  as  interparticle  sliding,  resulting  in  a  greatly 
decreased  overall  shear  resistance  and  a  lower  strain  at  peak  resistance.  As  well,  systems  composed  of  more 
rounded  particles  exhibit  fabric  differences  such  as  a  decreased  range  of  void  ratios  between  the  densest  and 
loosest  conditions,  lower  average  void  ratios,  and  less  variation  in  void  ratio  throughout  a  sample. 

Micromechanically-based  constitutive  models  based  on  disks  and  spheres  which  explicitly  account  for 
particle  rotaticm  have  recently  been  developed.  Disk-based  and  sphere-based  Discrete  Element  Method  (DEM) 
computational  models  are  common  in  micromechanical  and  geotechnical  modeling.  Since  these  models  may 
exhibit  considerably  more  rolling  behavior  than  is  found  in  real  (ie  more  angular)  soils,  it  is  particularly 
impotant  to  quantify  the  effect  of  particle  shape  and  the  rolling  mechanism  on  granular  material  behavior  at  this 
point  in  time. 

This  Repot  presents  the  results  of  an  Air  Force  Office  of  Scientific  Research-funded  study  on  the  effect  of 
particle  sh^  on  the  overall  material  behavior  of  granular  materials  through  a  Discrete  Element  numerical  study 
using  ellipse-shaped  particles.  The  Discrete  Element  Method  possesses  unique  capabilities  for  studying  the 
behavior  of  particulate  materials.  Because  it  is  a  numerical  simulation,  internal  contact  forces,  particle 
cxientatiois  and  average  stresses  and  strains  may  be  monitored  inobstrusively,  and  reproducible  "samples"  can 
be  guaranteed  from  test  to  test.  Because  the  DEM  models  each  particle  explicitly,  particle-to-particle  contacts 
nu^  be  bix^coi  and  roiewed,  allowing  fcx  large  local  deformations.  This  stu<fy  involved  calibrating  a  previously 
developed  two-dimensional  ellipse-based  computational  model  (Ting  and  Corkum,  1993)  with  physical  test  data 
by  Oda  et  al  (1983)  on  oval  cross-section  rods  and  by  Chapuis  (1976)  on  circular  rods. 

Additional  numerical  simulations  were  then  conducted  to  determine  the  relative  importance  of  angularity, 
contact  friction,  form  and  magnitude  of  contact  stiffiiess,  particle  gradation  and  degree  of  densification. 
Numoical  sanq)Ies  ccxisisting  of  multiple  sized  particles  of  varying  particle  aspect  ratio  were  formed  by  packing 
particles  with  their  major  axes  along  preferred  bedding  orientations  ranging  from  horizontal  to  vertical  bedding. 
Individual  particle  aspect  ratio  (major  axisminor  axis)  ranged  from  1:1  to  3:1.  The  samples  were  initially 
compressed  to  isotropic  conditions,  then  sheared  in  biaxial  shear.  This  Report  presents  the  results  from  these 
tests  along  with  the  conclusions  which  were  drawn  based  on  these  data. 
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2.  BACKGROUND 


One  of  the  most  impcxtant  facbxs  governing  the  behavior  of  cdiesirailess  granular  materials  is  particle  shape. 
This  Sectiai  details  the  effects  of  particle  shape  on  strength,  densification,  pore  pressure  buildup  and  mechanisms 
ofdeformaticHi. 

The  mechanical  behavior  of  soils  is  a  function  of  soil  structure  and  applied  stress  system.  Soil  structure  is 
ccxnposed  of  two  axiqxmoits  fabric  (spatial  arrangement  of  particles)  and  nature  and  magnitude  of  interparticle 
fcxces.  Fw  cdiesirmless  soils,  the  interparticle  forces  are  largely  governed  by  the  fabric  (Ladd  et  al  1977).  The 
four  factors  whidi  are  generally  accept^  as  having  the  most  influence  on  the  fabric  (and  hence  the  behavior)  of 
a  cohesionless  particulate  material  are:  particle  shape,  relative  density  D, ,  gradation  and  mineralogy. 

Of  these  factors,  particle  shape  and  relative  density  are  particularly  important.  Terzaghi  and  Peck  (1967) 
report  an  increase  in<))of5toll°fcH’  sands  with  increasing  angularity,  while  Lambe  and  Whitman  ( 1 968)  report 
a  5  to  6°  increase.  Holubec  and  D'Appolonia  (1973)  reported  an  increase  in  <|)  of  10°  by  increasing  angularity 
(see  Figure  1)  and  5°  by  increasing  £),..  It  is  often  difficult  to  separate  the  effects  of  these  two  factors,  as  the 
particle  shape  affects  the  range  of  void  ratios  and  possible  for  a  given  gradation,  and  hence  influences 
the  actual  density  range  over  which  a  particular  relative  density  is  computed.  Koemer  (1970)  found  that 
mineralogy  was  impcxtant  to  the  observed  internal  fiicticsi  angle  <)>  throug)i  its  influence  on  the  physical  properties 
of  the  particles,  ie  shape  and  interparticle  fiiction  Several  researchers  report  increases  in  4>  as  the  material 
becomes  less  unifcxm  in  particle  size  {ie  becomes  well-graded).  Lambe  and  Whitman  (1968)  report  increases  of 
2  to  4°  between  uniform  and  well-graded  materials. 

Particle  sh^  is  often  quantified  by  the  terms  sphericity  and  roundness  (Feda  1982),  where  sphericity  refers 
to  how  close  the  overall  particle  shape  is  to  a  sphere  and  roundness  refers  to  how  rounded  the  comers  are.  The 
coefiQcient  of  angularity  E  used  in  Figure  1  is  a  function  of  the  reciprocal  of  the  sphericity,  with  a  value  of  umty 
for  spheres  and  increasing  with  increasing  particle  angularity. 

2.1  Effect  ofParticle  Shape  on  Strength 

The  results  in  Table  1  indicate  that  increasing  particle  angularity  can  increase  measured  material  angle  of 
internal  firicticm(|)  by  about  5°.  Terzaghi  and  Peck  (1967)  report  increases  in  (|>  of  up  to  11°.  Koemer(1970) 
reported  increases  in  <|>  of  6  to  8°  with  increasing  angularity.  Holubec  and  D'Appolonia  (1973)  ran  tests  on 
matflrials  with  similar  coefficients  of  uniformity  and  coefficient  of  curvature .  Figure  1  plots  (j)  as  a  function  of 
angularity,  and  shows  that  4»  varies  up  to  10°  over  the  range  of  angularities  tested,  where  the  "angularity”  was 
computed  fi'om  specific  surface  areas  determined  indirectly  fi’om  permeability  tests.  In  addition,  Holubec  and 
D'Appolonia  report  that  the  axial  strain  at  peak  strength  €p^  increased  with  increasing  angularity.  For  glass 
beads,  almost  no  change  was  observed  in  while  the  angular  sands  show  considerable  variation  in  €p^  with 

D. 

2.2  Effect  ofParticle  Shape  on  Densification 

According  to  Youd  (1973),  the  primary  factors  controlling  the  maximum  and  minimum  void  ratio  limits  of 
sands  are  particle  shape,  particle  size  range  and  the  shape  of  the  gradation  curve.  His  data  indicate  that 
increasing  roundness  leads  to  increasing  densities  (lower  void  ratios)  as  well  as  a  decreased  range  of  density  - 

ThesetrendswCTec(MifirmedbyHolubecandD'Appolonia(1973) andDickin(1973).  Bhatia and Soliman 
(1990)  analyzed  the  distribution  of  void  ratios  fi’csn  secticxis  of  samples  and  found  that  the  range  in  the  local  void 
ratios  increased  {ie  greater  local  density  variation)  as  the  particle  shape  became  more  angular  and  the  relative 
density  increased. 
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Coefficient  of  Angularity,  E 

Figure  1 .  Effect  of  particle  angularity  on  peak  angle  of  internal  friction  (fi'om  Holubec  & 
d’Appolonia  1973). 

23  Effect  of  Particle  Shape  on  Pore-pressure  Buildup 

Accoding  to  Ishibashi  etal  (1982),  the  soil  properties  having  the  greatest  influence  on  liquefaction  are;  (1) 
grain-size  distributicHi,  (2)  soil  imifotinity,  (3)  grain  shape,  and  (4)  void  ratio.  According  to  them,  properties  (1) 
and  (2)  may  be  represented  by  the  mean  particle  diameter  and  coefficient  of  imiformity  C„ ,  (3)  by  the 
sphericity  airf  (4)  by  the  natural  void  ratio  e  and  Their  data  indicate  that  the  applied  shear  stress  ratio  t/o,,' 
to  cause  liquefaction  at  the  30th  cycle; 

•  increased  with  increasing  angularity 

•  increased  with  increasing  C,  up  to  C=5 

•  increased  with  decreasing  volume  change  potential  e  - 

In  other  wixds,  increasing  angularity,  decreasing  sample  uniformity  and  decreasing  the  potential  volume  change 
all  tend  to  increase  the  resistance  of  the  material  to  potential  liquefaction. 

These  results  woe  cxaiflrmed  by  Vaid  et  al  (1985),  who  also  found  that  at  high  confining  stresses,  samples 
madfi  of  angular  particles  can  still  be  made  to  Uquefy,  ^iiile  samples  made  of  rounded  particles  tended  to  stabilize 
above  a  relative  density  of  85%.  Thty  attributed  this  to  contact  crushing  and  breakage  of  the  angular  particles 
at  the  high  confining  stresses. 

2.4  Effect  of  Particle  Shape  on  Mechanisms  of  Strength  and  Deformation 

Based  on  the  deformational  behavior  of  regularly  packed  rigid  spheres,  Rowe  (1962)  proposed  his  stress- 
dilatancy  theory  wherein  he  decomposed  the  shear  resistance  of  cohesionless  granular  materials  into  three 
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components: 

•  fhctional  resistance,  due  to  interparticle  sliding, 

•  interference  resistance,  where  particles  rearrange  and  reorient  at  constant  voliune,  and 

•  dilatant  resistance,  due  to  the  effect  of  dilation  against  a  confining  stress. 

Rowe's  stress-dilatancy  theoiy,  and  Hcxne's  (1965, 1969)  revisions  did  not  consider  the  possibility  of  rolling 
between  particles.  Skinner  (1969)  showed  that  for  an  assemblage  of  uniform  spheres,  the  rolling  phenomenon 
is  voy  impotant.  By  direct  shear  tests  (xi  diy  and  wet  glass  ballotini,  steel  balls  and  lead  shot,  Skinner  measured 
peak  <|)„„  values  varying  firan  30  to  40°,  depending  on  initial  porosity.  However,  the  constant  volume  angle  (j)„ 
was  observed  at  about  25  °  regardless  of  <|)^,  even  though  <1>^  varied  from  4°  to  40°.  Skiimer  concluded  that  as 
the  sliding  fiiction  angle  increased,  the  sliding  mechanism  became  less  prominent,  and  rolling  predominates. 

This  conclusion  is  confirmed  by  tests  on  two  dimensional  photoelastic  oval  rods  (Oda  et  al  1983,  1985, 
Konishi  et  al  1983,  Mehrabadi  et  al  1988).  These  were  biaxial  compression  loading  tests  which  included  two 
diffoent  particle  aspect  ratios  (1.1  and  1.4),  two  interparticle  fiiction  angles  <|)^  (26°  and  52°)anda  variety  of 
bedding  angles  (0°,  30°,  60°  and  90°).  By  considering  the  statistics  of  the  contacts  where  pure  rolling,  pure 
sliding  and  both  rolling  and  sliding  were  occuiring,  as  well  as  the  contacts  where  the  fiictional  resistance  was  fully 
mobilized,  they  concluded  that  the  rolling  mechanism  dominates  when  the  interparticle  fiiction  angle  is  large. 

In  particular,  note  the  samples  of  nearly  round  particles  (aspect  ratio  1.1,  Oval  1)  exhibit  relatively  small 
dependence  on  (j),,  or  bedding  plane.  For  the  more  angular  particles  (aspect  ratio  1.4,  Oval  11),  there  is  a 
pronounced  drop  in  shear  strength  at  the  bedding  angle  of  60°,  as  expected  when  the  plane  of  maximiun  t/o 
closely  aligns  with  the  bedding  plane.  However,  since  the  material  with  the  higher  <)>^  exhibits  a  much  higher 
strength,  this  suggests  that  sliding  is  more  important  for  the  more  angular  particles.  Conversely,  this  also 
suggests  that  rotation  is  inhibited  for  the  more  angular  particles.  Direct  comparison  of  the  strengths  for  Oval  1 
and  Oval  n  at  the  same  fiiction  angle  (t>^  at  the  0°  bedding  angle  indicates  that  the  more  angular  particles  exhibit 
more  than  double  the  strength  of  the  more  rounded  particles. 

The  conclusicxi  that  particle  rolling  is  an  impcxtant  deformational  mechanism  is  further  supported  by  various 
experimental  and  numerical  data  on  two  dimensional  circular  and  oval  rod  systems.  Chapuis  (1976)  tested 
sanqjles  of  circular  rods  and  found  (J)„  of  15  °  and  <|),„„  of  25  °  for  a  dense  sample.  They  measured  the  rod-surface 
fiicticxi  angle  <j>^  as  15  to  16°  (greater  than  their  <|>„)  based  on  shearing  parallel  to  the  axis  of  the  rods  {ie  sliding 
only)  and  10  to  1 1  °  when  shearing  occuned  perpendicular  to  the  direction  of  the  rods  {ie  rolling  and  sliding). 

Ting  et  al  (1987)  simulated  the  shearing  of  two  dimensicmal  disk  systems  using  the  Discrete  Element  Method 
(DEM).  The  DEM  is  a  numerical  method  which  treats  each  particle  distinctly,  allowing  interparticle  slip  and 
rolling,  and  will  be  described  in  greater  detail  later  in  this  Proposal.  They  reported  the  results  of  shearing  disks 
indirect  sin:q)le shear.  For  an  interparticle  fiiction  <j>^  of  0°,  the  observed  ((>„„  was  5°,  while  for  =  25°, 
was  19°.  In  simulations  of  "triaxial"  tests  with  rigid  and  flexible  side  boundaries,  was  20°  for  4>^  =  25°, 
correlating  well  with  the  results  fi'om  the  simple  shear  simulations.  The  individual  particle  rotations  could  be 
inhibited  by  artificially  increasing  the  individual  particle  moments  of  inertia  I^.  When  no  rotation  is  permitted 
{ie  pure  translation),  the  peak  observed  firiction  angle  could  be  increased  up  to  40°  for  <|)j,  =  25  °. 

All  the  data  cited  in  the  previous  section  suggest  that  rolling  can  be  an  important  deformational  mechanism 
in  granular  materials,  especially  wiien  the  particles  are  round.  However,  real  soils  are  far  from  the  nearly  perfect 
sphoes  (x*  disks  considered  in  most  of  these  tests.  As  a  result,  real  soils  may  actually  exhibit  much  less  rolling 
that  might  be  suggested  from  the  above  data.  In  the  stucfy  by  Oda  et  al,  the  aspect  ratio  of  the  ovals  was  varied 
only  from  1.1  to  1.4  (equivalent  to  roundness  values  of  0.91  to  0.71  and  sphericity  values  of  0.998  to  0.97, 
respectively).  Nevertheless,  they  observed  a  twofold  increase  in  peak  resistance  with  decreasing  roundness. 
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3.  RESEARCH  PROGRAM 


The  data  cited  in  the  preceding  Section  has  shown  that  particle  shape  critically  influences  the  stress-strain- 
strength  behavior  of  drained  and  undrained  cohesionless  particulate  materials.  It  is  hypothesized  that  the  much 
of  the  observed  behavior  differences  between  rounded  and  angular  materials  is  due  to  a  different  packing  and 
densification  behavior,  and  altered  controlling  deformational  mechanisms.  In  particular,  samples  composed  of 
rounded  particles  more  readily  permit  interparticle  rolling,  with  a  resultant  weaker  observed  shear  resistance. 

In  the  research  areas  relating  to  particulate  mechanics,  various  analytical  and  numerical  models  based  on 
disks  and  spheres  have  emerged  and  become  fairly  popular.  Recently,  Chang  and  Liao  (1990)  developed  an 
analytical  micromechanically-based  constitutive  relationship  for  disks  and  spheres  accounting  for  partible 
rotaticHi.  These  ccmtinuum  formulaticxis  may  be  incorporated  into  finite  element  code  to  incrementally  solve  full 
boundary  value  problems.  This  voy  premising  development  must  be  tempered  by  the  knowledge  that  "real"  soils 
possess  hetoo^nous  strain  fields  (Chang  and  Misra  1990)  and  also  may  actually  roll  a  great  deal  less  than  that 
predicted  by  any  round  particle  model. 

Numerically,  the  Discrete  Element  Method  (DEM)  models  particulates  as  a  collection  of  individual  particles, 
and  was  originally  developed  by  Cundall  (1974)  to  model  rock  masses  composed  of  angular  blocks.  This  has 
spawned  various  two-dimaisicmal  disk-based  programs  such  as  BALL  (Cundall  and  Strack  1979),  DISC  (Corkum 
and  Ting  1986),  D/SC  (Bathurst  and  Rothenburg  1988),  SKRUBAL,  CONBAL2  (Dobry  and  Ng  1989)  and  three- 
dimensional  sphere-based  programs  TRUBAL  (Cundall  1988),  3DSHEAR  (Walton  et  al  1988),  DMC  (Taylor 
and  Preece  1 989),  and  SPHERE  (Hassan  1 990).  Other  disk-bas^  algorithms  have  also  been  developed  (Kishino 
1988,  Bazant  et  al  1990).  These  numerical  simulations  appear  to  quantitatively  approximate  the  behavior  of 
disks  and  spheres  appropriately  (Ishibashi  et  al  1989,  Chen  1989).  However,  approximation  of  realistic  soil 
behavior  using  disk-bas^  models  (especially  at  fairly  large  strains)  requires  inhibition  of  individual  particle 
rotation  (Ting  et  al  1989). 

Because  of  the  premise  exhibited  by  the  analytical  as  well  as  numerical  developments  in  particle  modeling, 
it  was  considered  to  be  especially  important  to  investigate  the  effect  that  using  a  circular  particle  form  factor 
would  have  <»  the  ovwall  response  of  a  particulate  system.  In  particular,  to  what  degree  does  out-of-roundness 
affect  material  behavior?  To  gain  insight  to  these  and  other  questions,  a  previously-developed  ellipse-based 

implftmmtatinn  of  the  Disoete  Elonoit  Method  was  used.  The  DEM  possesses  unique  capabilities  for  stuefying 
the  bdiaviw  of  particulates,  firm  its  ability  to  model  the  fine  detail  of  mechanisms  of  deformation  to  generating 
the  face  and  orioitation  distributions  during  numerical  testing.  The  computer  model  was  first  used  to  simulate 
the  physical  experiments  of  Oda,  Konishi  and  Nemat-Nasser  on  oval  rods  and  Chapuis  (1976)  on  circular  rods 
by  varying  the  contact  parameters  and  contact  stiffhess. 

For  the  second  stage  of  the  research,  numerical  samples  of  approximately  700  ellipse-shaped  rods  were 
formed  with  the  particles  aligned  along  prefened  bedding  planes,  then  isotropically  compressed,  and  sheared  in 
biaxial  compression.  Samples  were  tested  with  bedding  planes  varying  from  0,  30, 60  and  90°  normal  to  the 
maja  principal  stress  directiai,  as  well  as  "random"  orientation  (no  preferred  bedding).  Different  samples  were 
famed  at  varying  particle  aspect  ratio,  with  the  ratio  of  the  majoriminor  axes  of  the  individual  particles  varying 
fron  1 : 1  to  3: 1,  with  the  aspect  ratio  held  constant  in  each  sample.  As  well,  interparticle  angle  of  sliding  fnction 
(J)^  sample  size  (number  of  particles),  and  contact  prq)aties  woe  varied  in  some  instances.  This  Report  presents 
the  results  of  this  study. 
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4.  DISCRETE  ELEMENT  METHOD 


The  Discrete  (or  Distinct)  Elemoit  Method  (DEM)  has  been  used  for  over  a  decade  for  numerically  modeling 
the  mechankal  bdiavior  of  granular  materials.  Instead  of  a  continuum,  the  DEM  treats  granular  material  as  an 
assemblage  of  distinct  particles,  each  governed  by  physical  laws.  Each  particle  interacts  with  its  neighbors 
through  particle-to-particle  contacts  which  can  be  formed  or  broken  at  each  time  step,  as  shown  in  Figure  2. 
Because  no  restrictions  exist  regarding  mode  of  deformation  or  amount  of  displacement  of  each  particle,  the 
DEM  is  uniquely  suited  for  modeling  large  defecation  processes  and  bifijrcation-type  behavior  in  granular 
materials. 

This  numerical  technique  was  fu^t  adapted  by  Cundall  (1974)  to  geomechanics  for  stuefying  the  dynamic 
behavior  of  rock  masses.  As  this  technique  has  evolved,  it  has  been  used  in  a  wide  variety  of  research 
applications  in  engineering  mechanics  and  geotechnical  engineering.  These  applications  include  statistical 
micmnedianics  (Cundall  and  Strack  1979,  Cundall  1988,  Bathurst  and  Rothenburg  1988),  constitutive  behavior 
of  granular  sc^  (Zhang  and  Cundall  1986,  Ting  and  Cotkum  1988,  Dobry  and  Ng  1989,  Issa  and  Nelson  1989), 
creq)  of  soils  (Kuhn  and  Mitchell  1989),  flow  of  granular  materials  (Walton  et  al  1988,  Campbell  and  Bronwen 
1983,  CHiaboiKsi  and  Barbosa  1988,  Hakuno  et  al  1980),  analysis  of  rock-support  interaction  (Lorig  and  Brady 
1984),  blast  loading  of  ground  (Butkovich  et  al  1988,  Taylor  and  Preece  1989)  and  large  deformation 
geotechnical  modeling  (Barbosa  and  Ghaboussi  1987,  Ting  et  al  1989). 

The  largest  problems  tackled  to  date  with  the  DEM  involve  up  to  a  few  tens  of  thousands  of  two-  or  three- 
dimensional  particles  simulated  with  a  few  tens  of  thousands  of  time  steps.  The  particles  are  usually  of  simple 
germietric  sh^,  either  two-dimoisional  discs,  polygons,  three  dimensional  spheres  or  polygonal  blocks.  While 
circular  particles  have  the  great  advantage  of  computational  simplicity,  thq^  also  possess  the  inherent  tendency 
to  roll,  h  particular,  ccxnputed  ncxmal  contact  forces  never  cemtribute  to  the  moment  acting  on  a  particle  as  each 
ncxmal  contact  force  always  acts  through  the  particle  centroid.  As  a  result,  it  is  difficult  to  simulate  "real" 
material  behavior  when  the  basic  particle  is  even  slightly  angular  (Ting  and  Corkum  1988). 


For  particle  i, 
contact  forces 
are  summed 
over  neighbors 
j,  k,  I,  m  and  n 


Figure  2.  Lcxsal  equilibrium  scheme  in  Discrete  Element  Method 

An  alternative  is  to  use  polygonal  particles,  which  can  present  difficulties  at  the  vertices.  For  instance,  one 
must  distinguish  between  edge-edge,  edge-comer  and  cmner-comer  contacts,  and  also  maintain  a  smooth 
transition  in  bdiavior  from  one  contact  type  to  another.  In  addition,  polygons  have  no  unique  outward  normal 
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at  each  vertex,  and  are  c(xnplex  to  extend  to  three  dimensional  blocks.  In  particular,  approximating  the  shape 
of  rounded  or  giihangiilar  granular  matoials  such  as  cohesionless  soil  or  powders  using  polygonal  particles  may 
require  using  a  large  number  of  sides,  hence  increasing  the  number  of  corners  and  sides  in  a  particle. 

Because  of  these  difficulties,  it  was  decided  to  implement  the  DEM  using  two-dimensional  ellipse-shaped 
particles.  The  ellipse  shape  has  the  advantage  of  having  a  unique  outward  normal  and  no  singularities  at  every 
pmntalcxig  its  surface.  Extoisiai  to  dueedimoisiaisdiould  be  relatively  simple.  The  ellipse  falls  into  a  general 
class  of  fiincticxis  known  as  siqierquadrics  have  beoi  prc^xised  for  modeling  solids  (Williams  and  Pentland 

1989).  However,  because  the  ellipse  has  a  maximum  degree  of  only  two,  roundoff  errors  are  minimized  in 
compaiiscxi  with  s(xne  of  the  higher  degree  supoquadnc  fiincticMis.  Because  of  these  properties,  and  the  inherent 
limitations  of  using  circular  particles  for  modeling  "real"  granular  materials,  it  has  recently  become  popular  to 
use  the  ellipse  shape  for  Discrete  Element  modeling  (Tingl991,  Wei  ero/ 1991,  Rothenberg  and  Bathurst  1991, 

Ng  1992). 

This  Sectim  presents  the  detailed  development  of  the  numerical  algorithm  for  the  Discrete  Element  Method 
using  two-dimensional  ellipse-shaped  particles.  This  Section  derives  the  equations  for  determining  particle- 
particle  and  particle-wall  contact  locations  and  incremental  velocities  and  integrates  them  into  the  conventional 
DEM  algorithm.  Initial  validation  tests  on  rod  assemblies  of  varying  aspect  ratio  are  presented. 


Figure  3.  Nature  of  particle-particle  contact  in  Discrete  Element 
Method. 


4.1  Derivation  of  Ellipse-based  Discrete  Element  Method  Algorithm 

In  the  DEM,  every  particle  in  the  soil  mass  is  identified  separately,  with  its  own  mass,  moment  of  inertia  and 
contact  prcqierties,  as  shown  in  Figure  2.  Each  particle  (or  discrete  element)  is  considered  to  be  geometrically 
inviolate,  with  deformable  contacts.  Normal  and  tangential  springs  and  dashpots  exist  at  each  contact  to 
represent  ccmtact  compliance  and  energy  absorption,  with  a  limit  placed  on  the  maximum  tangential  shear 
strength,  as  seen  in  Figure  3.  The  contact  springs  may  be  linear  or  nonlinear,  and  in  some  cases  have  been  based 
oa  Mindlin’s  solutions  fix  Hotzian  contact  (Cundall  1988,  Dobiy  and  Ng  1989).  The  amount  and  rate  of  overlap 
between  neighboring  elements  are  used  to  determine  the  contact  forces  at  each  instant  in  time.  The  total 
unbalanced  fixces  and  moment  acting  on  each  particle  are  computed  based  on  local  equilibrium  of  direct 
neighbors  and  are  used  to  estimate  each  particle's  accelerati<Mis.  The  accelerations  are  then  integrated  for 
velocities  and  displacements  using  a  second  order  leapfrog  scheme  (Cundall  and  Strack  1979). 
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for  all  time  do 

j  store  data  to  fHes,  If  necessary 

j  conipiite  flexible  wall  stresses,  necessary 

I  forallpartidesdo 

j  I  for  all  items  in  contact  with  current  particle  do 

I  I  j  compute  contact  location  and  veloc¬ 

ities 

I  I  I  compute  contact  forces 

j  I  j  update  particle  and  veill  forces 

I  I  compute  accelerations  for  current  particle 

j  update  rigid  wait  velocities 

I  for  all  particles  do 

j  I  integrate  for  velocities  and  displacements 

j  j  form  list  of  particles  to  check  for  new  contacts 

j  :  for  all  walls  do 
j  I  update  wall  locations 


TABLE  1.  Pseudocode  for  Discrete  Element  simulation. 


The  flow  diagram  fw  the  g^eral  DEM  algorithm  is  shown  in  pseudocode  in  Table  1.  The  following  Sections 
present  the  derivation  of  the  ellipse-based  DEM  algorithm  based  on  computations  of  ellipse-ellipse  and  ellipse- 
wall  intersection,  and  associated  incremental  contact  velocities  (Ting  1991). 


4.1.1  Ellipse-ellipse  Intersection 


Since  the  DEM  cxxisiders  particles  to  be  quasi-rigid  with  deformable  contacts,  one  can  arbitrarily  defme  the 
"point"  of  contact  (x^ ,  y^)  between  two  ellipses  as  the  midpoint  of  the  line  connecting  the  two  ellipse-ellipse 
intersecticxi  points.  To  solve  for  these  intersections,  note  that  the  equation  of  a  two-dimensional  ellipse  located 
at  the  origin  of  its  local  coordinate  system  is  simply; 


-1=0 


(1) 


where  a  and  b  are  the  major  and  minor  axes  aligned  along  the  x  andy  coordinate  axes,  respectively. 
Transformation  between  local  and  world  coordinates  can  be  performed  using  the  following; 


JC 

cosd  -sin^  dx 

JC 

*  zz 

sin^  cos^  dy 

1 

\vorid 

0  0  1 

1 

(2) 


v^hCTe  ^is  the  orientaticm  of  the  local  system  with  respect  to  the  world,  measured  counterclockwise  positive.  For 
an  ellipse  with  center  at  fcbc,  dy)  in  world  coordinates  and  major  axis  rotated  at  angle  5  to  the  x  axis,  measured 
counterclockwise,  the  general  equation  in  world  coordinates  is; 

fix,y)  =  A(x-dxf*B(y-dy)^^2C(x-dx)(y-dy)-l  =  0  (3) 


where  coefficients  A,  B  and  C  are  functions  of  dx,  dy,  a,  b  and  ^only; 
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(4) 


A  = 
B  = 


'  COS0^ 

2 

'  sin^^ 

b  j 

(  sin^^ 

2 

^  cos^'' 

i  0  J 

^  b  J 

C  =  COS0  sin^ 


J_  _  _1_ 


Fot  two  ellipses /andy  with  axes  o,-,  A,  anda^,  bj,  centers  iXci,y„)  and  {x^  ,yg)  with  rotations  0^  and  Oj,  the 
^latinn  fcM"  ellipse  J  expressed  in  the  local  coordinate  system  of  ellipse  /  is  simply  equation  (3),  except  that  Oj, 
bj  replace  a  and  b  in  equation  (4),  and  the  new  6,  dx  and  dy  are  given  by: 

e=0j-0, 

^  cosl^.  +  iy^j  -  yj  sin^.  (5) 

^  =  -(^cj  -  +  (y^j  -  yJ  cos0, 


Since  particles  in  a  DEM  simulation  are  usually  only  marginally  in  contact,  it  is  important  to  minimize  the 
pfiffrf  of  aity  roundoff  etiw  in  ccxiputing  ellipse-to-ellipse  intersecdrais.  As  well,  the  intersection  equation  should 
maintain  the  same  accuracy  regardless  of  the  units  used,  and  should  allow  for  a  wide  range  in  the  magmtude  of 
the  aspect  ratio  a/b.  To  address  these  issues,  first  transform  the  ellipses  to  the  local  coordinate  system  of  ellipse 
I  using  equations  (3),  (4)  and  (5),  then  fiirther  transform  the  system  by  scaling  both  local  x  and  y  axes  to 
homogeneous  form; 


2 


(6) 


Figure  4.  Nomenclature  fw  contact  between  two  ellipses. 

Ellipse  /  then  reduces  to  a  unit  circle; 

-  I  =  0  (7) 

The  equation  for  ellipse  J  in  the  new  homogeneous  local  coordinates  X  and  T  becomes: 
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(8) 


A  (X  -  dxf  +  B  {Y  -  dyf  ^  2  C  {X  -  dx)  (Y  -  dy) 
where 

A  =  afA,  B  =  bfB,  C  =  a.b^C,  ^  ^  ^  (9) 

fl/  bf 

To  solve  for  the  intersection,  Y  is  eliminated  from  equations  (7)  and  (8),  resulting  in  a  quartic  equation  in 
toms  ofX  This  equation  has  beai  presented  elsewhere  (Ting  1991,  Rothenberg  and  Bathurst  1991).  However, 
it  turns  out  that  this  equation  can  be  quite  ill-condition^.  Since  the  particles  being  simulated  are  usually  fairly 
"hard",  the  amount  of  ovo*!^  is  usually  quite  small.  As  a  result,  the  location  of  contact  often  approaches  a  single 
point.  If  both  particles  are  fairly  flat  near  the  contact  location,  considerable  error  may  result  in  computing  the 
contact  location.  Also,  in  some  situations  where  both  particles  have  major  axes  aligned  with  each  other  (ie  6  -* 
0),  the  quartic  equation  nearly  reduces  to  a  quadratic,  Aat  is,  the  coefficients  oiX*  and  approach  zero.  This 
can  also  lead  to  (xoblems,  as  the  solution  process  requires  ncxmalizing  the  quartic  so  that  the  leading  coefficient 
of  A"*  is  unity.  If  this  leading  coefficient  is  nearly  zero,  then  the  other  coefficients  become  very  large  after  the 
normalizaticm,  resulting  in  considerable  roundoff  error. 

Alternatively,  it  is  better  to  solve  for  a  "contact"  location  for  the  two  ellipses  using  a  better-conditioned, 
quartic  polynomial  (Ting  1992).  As  seen  in  Figure  5,  the  fimctiai  g(x,y)  =  c  is  a  family  of  ellipses  sharing  origi^ 
orientation  and  aspect  ratio  with  varying  axis  magnitudes  a  Jo  depending  on  the  value  of  the  constant  c.  It  is 
possible  to  locate  a  function  g(x,y)  =  c  which  just  touches  the  ellipse  =  0  (at  point  A').  At  the  common 
point  between  these  ellipses /(5c, =  0  and  g(x,y)  =  c,  the  slope  of  each  ellipse  is  the  same.  Similarly,  <me  can 
locate  a  function  f(x,y)  -  c'  wUch  just  touches  ellipse  g(x,y)  =  0  (at  point  B').  It  is  possible  to  define  a  new 
function  h(x,y)  =  0  which  is  the  locus  of  all  points  where  f(x,y)  =  c  and  g(x,y)  =  c'  share  the  same  slope.  This 
function  h(x,y)  =  0,  shown  in  Figure  5,  intersects  both  f(x,y)  =  0  and  g(x,y)  =  0  nearly  at  right  angles.  The 
solution  of  the  h(x,y)  =  0  with  g(x,y)  =  0  and  h(x,y)  =  0  vn(lif(x,y)  =  0  results  in  two  quartic  equations  which 
are  better  conditicxi^  than  the  quartic  resulting  from  merely  intersecting  two  ellipses.  The  resulting  intersections 
can  then  be  combined  to  locate  a  "point"  of  contact,  just  as  with  ellipse-ellipse  intersections,  only  now  the 
presumed  contact  location  (x^  ,yc)  is  midway  between  points  A'  and  B'. 


ellipse  j 
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Locating  the  intersection  of  h(x,y)  =  0  with  i  (x,y)  =  0  (point  A*  in  Figure  5)  requires  the  solution  of  the 
following  quartic  in  terms  of  (derived  in  Ting  1992): 

{aP  +  Q^)X*  +  (4PiR_+  2SQ)p  +  {P  -  aP  - 
{-2PR  -  2SQ)X  +  (P^  -  5^)  =  0 

w^ere 

P  =  C;  Q  =  B  -  A  ;  R  =  -B^  -  Cic ;  5  =  +  Crf 

Serving  equaticxi  (10)  results  in  iq)  to  fixir  roots  where  ■  (xy)  —  0  and  h(x,y)  =  0.  The  roots  are  tested  to  see  which 
has  the  mngnihiHe  ^(x,y)\  to  see  which  lies  closest  to  the  intersection  point.  The  local  coordinates  of  the 

intersection  arc  then  converted  to  world  cocH'dinates.  The  entire  system  is  the  solved  fOT  ellipse  j  in  terms  of 
ellipse  I  to  yield  a  second  intersection  which  would  correspond  to  Ae  intersection  of  g(x,y)  =  0  and  h(x,y)  =  0 
Ox)int  B').  The  presumed  locatiai  of  the  contact  is  then  taken  as  the  average  of  the  two  "intersection"  points  A' 
and  B',  that  is,  the  midpoint  of  the  line  joining  points  A’  and  B'.  The  normal  at  the  contact  is  also  taken  as  the 
average  of  the  outward  normals  to  each  ellipse  at  each  point  A'  and  B'. 

This  technique  results  in  a  much  more  reliable  and  accurate  algorithm  for  computing  ellipse-ellipse 
intCTsectiai  (Ting  1992).  In  additiai,  it  should  be  possible  to  extend  the  method  to  three  dimensions  more  easily 
than  computing  the  actual  surface  of  intersection  between  two  ellipsoids. 


(10) 

(11) 


4.1.2  Ellipse-Ellipse  Incremental  Contact  Velocities 


To  conq)ute  the  incronental  contact  fixee  betweoi  particles  I  and  j,  one  needs  to  estimate  the  relative  velocity 
of  the  contact  point  {x„  y^  on  particle  I  with  respect  to  the  same  point  <mi  particle/  Figure  6  shows  the  sign 
omventicMi  used  for  the  normal  and  tangential  vectors  at  the  contact,  where 


e  =  ej  +  ej 

A  A  A 

t  =  e^i  -  e^j 


(12) 


Figure  6.  Sign  convention  for  vectors  at  a  contact. 


For  particle  1,  the  magnitude  of  the  moment  arm  fiom  the  centroid  to  the  contact  is: 
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(13) 


in  local  or  global  coordinates,  while  the  radius  vector  in  the  direction  of  the  contact  is: 

*  (>'c  -  >’«•)/  =  '■c/COSfl'.f  +  r^^s\na.i 


The  relative  velocity  of  the  contact  point  on  I  relative  to  j  is  simply: 
dv  =  V .  -  V  . 

C  Cl  CJ 

where  the  velocity  of  the  contact  on  particle  I  is: 

=  (■'’a  -  ^ 

Then,  the  relative  velocity  of  the  contact  between  /  and  j  is: 

dv^  =  dvc^i  +  dvCyj 

=  K^x/  -  -  (v^j  -  d.r^jSina.)]!^  + 

-  (Vyj  +  d.r^jcosaj)/ 


(14) 

(15) 


(16) 


(17) 


This  relative  velocity  is  then  resolved  parallel  and  perpendicular  to  the  contact  normal  to  yield  the  incremental 
normal  and  tangential  contact  velocities: 

^'^cn  =  dvj  =  i{dv;  i)  =  (dvc^  +  dvc^  e^)  i 
^cs  =  dvj  =  tid\'i)  =  (dvc^e^  -  dvCye^)i 

Note  that  the  relative  contact  velocity  from  j  to  I  has  identical  magnitudes  dv„  and  dv,  as  equation  ( 1 8),  but  the 
outward  normal  and  tangential  vectors  differ. 


The  incremoital  contact  forces  during  a  given  time  step  may  then  be  computed  based  on  the  contact  spring 
stiffriess  k„  and  k„  contact  dashpot  values  c„  and  c,  and  time  step  of  integration  At: 


AF„  =  kdv^iit  +  c^dv^ 

ft  n  ft  ft  ft 

AF  =  kdvAt  +  cdv^ 

s  s  s  s  s 


(19) 


The  total  contact  force  may  be  computed  at  time  t  using  the  known  total  contact  force  at  time  t-At  using: 


-  AF. 

*  AF 

s 


(20) 


To  rq)resait  nonlinear  contact  behavior,  tangent  stiffiiesses  which  are  history-dependent  could  be  used  for 
the  normd  and  shear  spring  values.  The  maximum  tangential  shear  force  at  the  contact  is  usually  limited  by  a 
function  which  is  dependent  upon  the  normal  contact  force. 


It  turns  out  that  summing  incremental  contact  forces  at  each  time  step  can  lead  to  significant  accumulated 
OTOT.  This  was  noted  by  Ting  etal  (1989)  in  their  work  with  disk-shaped  particles.  Instead,  they  computed  total 
normal  overlap  at  each  instant,  then  used  this  total  overlap  to  obtain  the  total  contact  normal  force.  Tangential 
forces  were  still  obtained  using  incremental  forces,  as  in  Cundall  and  Strack's  original  algorithm.  For  ellipse- 
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ellipse  intersection  using  the  modified  intersection  routine  described  above,  it  is  possible  to  estimate  the  total 
normal  ovCTlap  as  the  distance  between  points  A'  and  B'  in  Figure  5.  These  points  are  merely  the  intersection 
oig(x,y)  =  0  with  h(x,y)  =  0  mdif(x,y)  =  0  with  h(x,y)  =  0,  respectively. 

This  procedure  is  carried  out  for  all  contacts  between  all  particles.  Then  the  contact  forces  are  summed  for 
all  contacts  k  for  each  particle  I  to  yield  net  particle  accelerations  at  time  t: 

=  E  [(•»,  -  -  ‘iFJ  *  (y,  - 


where  /w,  and  1  are  the  mass  and  polar  moment  of  inertia  for  particle  I,  respectively.  Body  forces  such  as 
gravitational  effects  may  be  added  in  the  force  sum  for  Fy  in  equation  (2 1). 


Using  the  leapfrog  integration  scheme  proposed  by  Cundall  and  Strack  velocities  of  each  particle  I  are 

computed  at  the  next  half  time  step  t+ViAt: 


^t+^AAt 

d 


The  displacements  and  rotation  of  each  particle  I  are  then  obtained  by: 

=  ^,-A/ 

yitAt  ~  yf-Ai 


Overall,  this  ejqjlicit  integraticxi  scheme  is  secaid-wdo-  accurate,  and  is  generally  regarded  as  the  best  overall 
fix  accuracy,  stability  and  efficiency  for  many-body  simulation  (Hockney  and  Eastwood  1981).  Global  viscous 
or  acceleration-proportional  damping  may  also  be  incorporated  using  the  techniques  outlined  by  Cundall  and 
Strack  (1979)  and  Zhang  and  Cundall  (1986). 

4.1.3  Ellipse-line  Intersection  and  contact  velocities 

The  calculation  of  ellipse-line  intersection  is  required  when  computing  the  interaction  between  an  ellipse 
particle  and  a  wall.  To  compute  the  intersection  points,  first  transform  the  equation  of  the  line  into  the  local 
r/wHinatftg  of  the  ellipse.  Then,  cme  may  solve  a  quadratic  equation  based  on  the  intersection  of  an  ellipse  1  and 
wall i  in  local  coordinates. 

Figure  7  shows  the  detailed  nomenclature  for  an  ellipse  /  and  wall  j  intersection.  For  a  wall;  with  endpoints 
in  local  coordinates  of  ellipse  /  equal  to  (x^ ,  y,)  and  (x^ ,  y^,  the  equation  of  the  wall  in  the  homogeneous  local 
coordinates  of  ellipse  /  is; 
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(24) 


y  =  m,X  *  b 


h  m 


^x] 


while  the  equation  of  the  ellipse  is  the  same  as  in  equation  (7).  Eliminating  Y  from  equations  (7)  and  (24)  results 
in  the  following  quadratic; 


(mf+l)X^  +  2mib,X  +  (bf  -  1)  =  0 


(25) 


The  roots  of  this  quadratic  are  substituted  back  into  the  local  ellipse  equation  (7)  to  obtain  the  full  coordinate 
pairs.  PCX'  a  pair  of  intosection  points,  the  contact  l(xation  may  be  taken  to  be  the  average  location,  which  is  the 
midpoint  (m  the  wall  between  the  intersections.  This  contact  location  is  then  transformed  back  into  wotM 
coordinates  (x„  using  equations  (7)  and  (2). 


Once  the  intersection  location  is  determined,  the  incremental  contact  velocities  need  to  be  calculated.  The 
wall  is  specified  by  a  ostAsx  of  lotatiai  (x^  and  endpoints  (x,,  y,)  and  (x^,  y^.  In  addition,  the  wall  possesses 
a  "sense"  which  is  ccmsidered  positive  when  the  solid  side  of  the  wall  is  found  clockwise  with  respect  to  the 
tangential  vector  Derived  quantities  include  the  angle  of  orientation  Uj  with  respect  the  horizontal,  and  the 
inM'nn/normalvectorf  and  tangential  vecto'il'from  endpoint  1  to  endpoint  2.  The  distance  from  the  wall  center 
of  rotation  to  the  contact  is  r„. 


(x2-y2) 


Figure  7.  Nomenclature  for  ellipse-wall 
intersection. 


The  wall  vectors  are  simply: 

H'  =  Wj  z*  +  Wj/ 

i  =  ej  +  ej  =  sense- {wj  -  wj ) 


(26) 
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The  velcxjity  of  the  ccxitact  cmi  ellipse  I  is  the  same  as  equation  (16),  while  the  velocity  of  the  contact  on  the  wall 


-  ('’J  *  V'>  ^  .  (27) 

=  ['W  -  ^j(y.  -  ^  [‘■j,  * 

where  and  are  the  components  of  the  velocity  vector  of  the  center  of  rotation  of  the  wall. 

The  relative  velocity  of  the  contact  from  ellipse  I  to  waliy  is  then: 
dv^  =  dvc^i  +  dvCyj 

=  [(V„  -  v^j)  -  6.(y^  -  yj  *  4(7,  -  7,,)]''  *  (28) 

[<v  -  >')»)  * 

This  relative  velocity  is  then  resolved  parallel  and  perpendicular  to  the  wall  in  the  same  manner  as  before,  with 
the  same  results  as  equation  (1 8). 

Resolving  the  normal  and  tangential  contact  forces  into  their  x  and  y  components  leaves  the  following 
contributions  to  the  force  and  moment  sums  for  wally  : 

F.  =  -  F^e, 

Fy  =  -F„e,  -  (29) 

=  F„[(x^  -  xje^  -  (y,  -  yje^] 

where  the  moment  is  taken  about  an  arbitrary  point  (x„  yj.  These  are  summed  for  all  particles 
adjacent  to  the  wall  to  obtain  the  total  wall  forces  and  moment. 

Use  of  the  inocmental  equations  (18)  for  wall-ellipse  contact  in  computing  the  total  contact  forces  can  also 
lead  to  significant  accumulated  error,  as  with  ellipse-ellipse  contact.  Alternatively,  the  total  overlap  may  be 
^proximated  as  the  ncxmal  distance  between  the  wall  and  the  point  on  the  ellipse  with  slope  parallel  to  the  wall 
(point  A  in  Figure  7).  Using  equaticm  (24)  as  the  equation  of  the  wall  in  local  homogeneous  coordinates  of  ellipse 
I  (for  m,  not  -  «=),  the  point  on  the  ellipse  of  common  slope  may  be  obtained  from: 

^  =  -L  =  -11  =  =  m  (30) 

dx  2Y  Y 

Substituting  (30)  into  the  equation  of  the  ellipse  (7)  and  eliminating  X  leads  to: 

Y  =  ±  1 - ;  X=  -mY  (31) 

^1  + 

After  converting  X,  Y  to  world  ccxjrdinates,  the  total  normal  overlap  4l/j  may  be  computed  by: 

^n=i>r^  (32) 

where  #  and  are  defined  in  Figure  8.  This  total  normal  overlap  is  then  used  directly  in  computing  total  normal 
contact  force. 
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4.2  Derivation  of  Associated  Properties 
4.2.1  Contact  properties 

The  typical  spring,  dashpot  and  slider  arrangement  at  each  contact  for  most  implementations  of  the  DEM 
is  shown  in  Figure  3.  In  the  curroit  fcsmulation,  this  has  beoi  modified  somewhat  in  that  it  is  presumed  that  each 
particle  possesses  a  set  of  springs  and  dashpots  at  each  contact,  which  combine  in  a  series  or  average  maimer. 
At  present,  normal  and  tangential  springs  and  viscous  dashpots  are  considered  to  act  in  series,  while  the  slider 
f(x  limiting  the  maximum  shear  strength  uses  an  average  set  of  strength  parameters.  This  allows  simulation  of 
variable  particle-particle  and  particle-wall  properties. 


The  spring  stiffiiesses  reccxnmended  for  use  can  be  inferred  fi'om  Hertzian  contact  theory  for  elastic  spheres 
oc  Q^inders  in  contact  (Cimdall  1988,  Ting  et  al  1989).  Alternatively,  one  may  use  empirical  modifications  to 
the  Hatz  values  to  account  fix-  local  yield  (Chang  et  al.  1989).  The  spring  stiffiiesses  recommended  for  use  can 
be  inferred  fixxn  Hertzian  contact  theory  for  elastic  spheres  or  cylinders  in  contact.  For  two  (^lindrical  disks  in 
contact  with  partial  slip,  Chang,  Misra  and  Xue  (1989)  quote  the  following  nonlinear  relationship  fw  normal 
spring  stiffiiess  k^. 


_ 2£G _ 

(1-t^  [2  lnClrlA)-\] 


(33) 


where 


'2r{\-v)F„Y 

nG 


(34) 


G,  V  are  the  shear  modulus  and  Poisson's  ratio,  respectively,  of  the  individual  particles,  r  is  a  combined  radius 
formed  from  the  individual  radii  r,  and  y  and  F  is  Ae  normal  contact  force.  In  the  current  computer 
implementatiiMi,  only  constant  spring  stiffiiesses  are  used. 


The  iKxmal  dashpots  at  eadi  oxitact  may  be  used  to  represent  the  energy  absorption  capability  of  the  particle. 
In  particular,  the  dasl^  value  nu^  be  directly  related  to  the  coefficient  of  restitution  e  for  the  individual  particle 
(Coikum  and  Ting  1986)..  The  combined  dashpot  value  c„  based  on  collision  of  two  particles  with  masses  m, 
and  m2 ,  with  combined  normal  spring  k„  is; 


The  maximum  shear  strength  at  the  contact  is  related  to  the  current  normal  spring  force  using  a  Coulomb-type 
fiiction,  as  indicated  in  Figure  3.  Tensile  strength  normal  to  the  contact  is  currently  set  to  zero,  although  efforts 
have  recently  been  made  to  allow  for  tensile  strength  in  an  effort  to  model  ice  interaction. 


vi^ere  m 


/n, 

w,  +/n2 


4.2.2  Stress  Tensor 

AvCTage  macroscc^ic  stresses  within  the  particle  assembly  may  be  estimated  using  the  external  wall  forces 
or  contact  forces  within  a  subregion  of  the  system.  Using  only  the  boundary  contact  forces,  the  average  stress 
tensor  in  a  volume  V may  be  computed  as  follows  (Cundall  and  Strack  1979): 
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1  C  T-C 
o..  =  —  V  Xi  Ti 


where  I,  j  are  direction  indices,  the  summation  occurs  for  all  contacts  c  on  the  boundary  S,  x,-  is  the  locaticm 
coordinate  and  Tj  is  the  boundary  traction  at  contact  c.  Alternatively,  the  average  internal  stress  tensor  m^ 
be  computed  using  all  the  contact  forces  within  the  specified  region  V  (Bathurst  and  Rothenburg  1988): 


where  the  summation  occurs  for  all  contacts  c  within  the  region  V,  p  is  the  contact  force  vector,  and  V  is  the 
contact  vector  (centroid-contact)  for  contact  c.  Note  that  for  circular  particles,  equation  (37)  degenerates  to 
equaticxi  (36)  as  all  CMitributiais  to  the  summatiai  from  interior  contacts  cancel  out  and  only  contributions  from 
the  surface  S  remain.  This  is  not  true  for  noncircular  particles,  however. 


4.2.3  Strain  Tensor 

Several  methods  have  been  used  for  computing  an  average  strain  tensor  from  a  deformed  assemblage  of 
granular  materials  These  methods  usually  rely  on  analyzing  the  deformation  of  an  external  membrane,  rather 
than  internal  displacements  of  the  particles.  Alternatively,  one  may  approximate  the  deformation  pattern  within 
a  region  with  a  constant  strain  finite  element  (Ting  and  Corkum  1992).  A  linear  displacement  field  may  be 
specified  by: 

V  =  pj  +  yfljJC.  +  P^y.  ^  ' 

\riiae  «,  V  are  the  displacements  in  the  x  and  y  directions,  respectively,  and  a,,  02,  0C2,Pj,fi2  are  constants 
to  be  determined.  The  "best"  strain  field  may  be  computed  from  particulate  data  by  using  the  best  least-squares 
fit  of  the  particulate  displacement  data.  The  sum  of  the  squares  of  the  error  for  each  particle  is: 

5,  =  E  Mx,,yp  -  uf  =  ^ 

(39) 

=  E  -  vf  =  EiA  Pih  ^ 

f=l  f=l 

where  yi)  are  the  initial  coordinates  of  particle  /,  14  and  are  the  x  and  y  displacements  for  the  particle, 
respectively,  and  the  summation  is  over  all  N  particles  in  the  region  of  interest.  Then,  taking  partial  derivatives 
of  equaticxi  (39)  with  respect  to  each  unknown  and  equating  to  zero,  the  following  two  systems  of  3x3  equations 
result: 

f  N  Ey,  fl,]  [  Ev,  1 


N 

Ea^, 

w 

Pi' 

E*, 

E^‘ 

P2 

= 

Ey, 

Ps, 

Solution  of  equation  (40)  for  the  variables  a,  ...Pi  results  in  the  following  for  the  average  strain  tensor  6j,: 
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43  Implementation  of  Ellipse-based  DEM 

The  ellipse-based  Discrete  Element  Method  algorithm  described  above  has  been  implemented  as  part  of  a 
menu-driven  intoactive  simulation  program  using  ANSI-standard  C  as  the  programming  language.  Many  of  the 
advanced  C(xicq)t5  fron  Cundall  and  Strack's  (1979)  BALL  program  are  implemented.  Single-linked  lists  were 
used  extensively  for  keeping  track  of  particles  in  contact  and  particles  in  close  proximity  to  one  another.  "Grid 
boxes"  or  "microcells"  are  used  to  assist  in  tracking  neighbcHS. 

At  every  instant  in  time,  contact  intersection  is  computed  between  every  particle  /  and  the  items  in  that 
particle's  "ccxitact  hst",  that  is,  those  balls  or  walls  which  are  in  contact  and  nearly  in  contact.  A  search  for  new 
ccxitacts  is  initiated  any  time  particle  /displaces  past  a  specified  threshold  value.  This  search  consists  of  checking 
for  contact  between  Ae  current  particle  and  the  current  contact  list,  as  well  as  the  particles  residing  in  the  grid 
box(es)  in  which  particle  I  resides.  This  contact  check  consists  of  successively  testing  for  overlap  of  the 
rectangular  boxes  circumscribing  each  ellipse,  then  overlap  of  the  circles  circumscribing  each  ellipse,  then  finally 
overlap  between  two  ellipses.  The  ellipse-ellipse  overlap  computation  is  performed  only  after  testing  positive 
on  the  simple  contact  checks.  Note  that  use  of  the  contact  and  grid  box  lists  reduces  the  order  of  the  problem 
from  an  "n-squared"  problem,  where  the  number  of  computations  is  proportional  to  the  number  of  particles  n 
squared. 

As  an  example.  Figure  8  plots  a  typical  system  of  particles  together  with  its  grid  box  scheme.  Particle  4 
(shown  in  dots)  resides  in  grid  boxes  6  and  7  and  is  drawn  together  with  a  boimding  ellipse  approximately  dn 
larger.  The  box  list  for  box  6  ccxitains  ellipses  1, 4, 5, 6  and  7,  while  box  7  contains  2, 3, 4  and  7.  While  particle 
4  is  actually  in  ccxitact  with  balls  1, 2, 5  and  6,  ball  7  also  exists  in  the  contact  list  for  particle  4  since  it  lies  within 
dn  of  ball  4,  even  though  these  particles  are  actually  not  in  contact.  At  every  time  step,  the  contact  between  ball 
4  and  balls  1,  2,  5,  6  as  well  as  7  are  recomputed  and  updated.  When  any  part  of  particle  4  moves  by  a 
"significant"  fraction  of  dn,  then  a  contact  check  is  initiated  for  all  the  particles  (and  walls)  lying  in  grid  boxes 
6  and  7,  ie  particles  1, 2, 3, 5, 6  and  7.  Note  that  wall  identification  numbers  can  also  be  included  in  the  grid  box 
list,  and  are  stored  as  negatives  to  distinguish  them  from  particles,  which  are  stored  positive. 
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Figure  8.  Contact  list  and  grid  box  list  schemes. 


Dynamic  memory  allocation  and  deallocation  occur  at  runtime  as  contacts  are  made  and  broken.  A  variety 
of  boundaries  are  si4}pofted,  including  rigid  stress-  and  face-controlled  walls  and  luiked  displacement-controlled 
Walk  Periodic  boundaries  are  being  in:q)lemented.  The  simulaticxi  program  is  text-based,  with  no  system-specific 
graphics,  and  has  been  ported  to  numerous  different  computers  such  as  IBM-ccanpatible  personal  computers, 
UNIX  workstatiois  and  minisupercomputers. 


4.4  Code  Validation 

The  was  initially  validated  using  several  simple  particulate  systems  to  ensure  compliance  with  physical 
laws  based  on  static  and  c^namic  equilibrium  coiditioiis  (Ting  and  Cakum  1991).  These  simple  test  simulations 
iiyhidod  particle-to-particle  and  particle-to-wall  collisions  as  well  as  systems  of  stacked  particles.  The  computed 
forces  and  energy  relationships  for  the  simple  test  simulations  compared  well  to  analytical  solutions. 

As  an  example,  Figure  9  plots  a  simple  close-packed  array  of  5  frictionless  disks  each  with  radius  1  m  and 
specific  gravity  2.67.  The  numbers  shown  on  the  left  half  of  the  Figure  are  theoretical  normal  contact  forces  in 
kN,  vriiile  the  italicized  numbers  on  the  right  half  are  computed  values  from  the  DEM  simulation.  The  normal 
Hamping  vvas  taken  as  the  value  which  would  result  in  a  coefficient  of  restitution  of  0.5,  with  disk-disk  and  disk- 
wall  nomal  spring  stiffiiesses  105  kN/m^  and  106  kN/n^,  respectively.  Timestcp  of  integration  was  0.00025 
sec.  The  sim^ilatinn  was  allowed  to  run  for  2  sec  (8000  timestepsX  and  d^amic  forces  were  observed  to  dissipate 
after  about  1  sec. 
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Figure  9..  Gn:q}aris(xi  of  theoietical  and  conputed  (italicized)  normal  contact  forces  for  a  sym¬ 
metric  firictionless,  close-packed  five  disk  system  subjected  to  gravity  loading.  Values  in  kN. 


4.5  Interface  Programs 

Also  implemented  are  separate  graphical  preprocessing  and  postprocessing  modules  incorporating  a  C  AD- 
s^e  interface.  The  graphical  interface  programs  allow  the  setup  of  granular  systems  for  simulation  and  allow 
the  viewing  of  ou^ut  files  written  the  simulaticxi.  These  personal  computer-based  programs  plot  force  chains, 
time  sequences,  animations,  contact  statistics  and  distributions,  and  also  compute  the  average  stress  and  strain 
according  to  the  algorithms  described  above.  These  are  described  in  greater  detail  in  Ting  and  Corkum  (1992). 

It  was  decided  to  maintain  sq^arate  (xograms  fcx*  the  simulation  and  the  interface.  This  allows  the  flexibility 
of  pcxting  the  main  con^ute-intensive  portion,  the  simulation,  to  the  fastest  computer  available  to  the  user.  The 
interface  programs  are  all  based  on  die  personal  computer-based  CAD  drawing  and  interface  package  FEINT 
developed  by  Corkum  (1990).  This  paradigm  of  separating  the  graphical  user  interface  (GUI)  from  the  main 
sdver  has  some  voy  real  advantages  not  (xily  in  its  current  application  with  the  DEM,  but  also  with  any  existing 
program  written  in  any  computer  language.  A  GUI  program  based  on  the  FEINT  package  may  be  tailored  to 
write  a  data  file  for  the  specific  solver  program.  Similarly,  a  tailored  post-processor  GUI  program  accepts  the 
output  file  fixxn  the  solver  program  f(x  graphical  interpretatitm  of  the  results.  This  removes  the  need  for  rewriting 
existing  batch-style  programs  while  making  available  the  easy-to-use  graphical  environment  and  mouse  support 
currently  enqjloyed  in  modem  coiiq)uter  programs.  At  present,  the  FEINT  graphical  interface  package  is  written 
in  ANSI-standWd  C  using  BorlandC-H-  and  also  Watcom  C-H-.  The  FEINT  package  consists  of  drawing 
functions  such  as  expected  of  a  two-dimensional  CAD  package,  such  as  zoom,  clip,  copy,  etc.  All  drawing 
entities  are  stored  throu^  calls  to  the  FEINT  library  in  a  drawing  database  maintained  using  a  double-linked  list. 
This  package  can  support  most  of  the  popular  IBM  PC-based  graphical  adapters,  with  automatic  hardware 
detection  available  at  runtime. 
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One  shortccHning  of  most  C-ianguage  packages  is  the  lack  of  a  plotting  capability  for  a  goierai  ellipse 
oriented  at  any  angle  to  the  horizontal.  Since  each  ellipse  particle  is  stored  in  the  drawing  database  using  the 
center  coo'dinates  and  angle  of  rotation,  and  since  the  FEINT  drawing  package  allows  zooming  into  parts  of  the 
drawing  database,  it  was  necessary  to  implement  a  pixel-based  plotting  routine  for  general  ellipses.  This  was 
done  by  modifying  Bresoiham’s  circle  algorithm,  where  decisions  on  plotting  are  made  on  a  pixel-by-pixel  basis 
(Foley  and  Van  Dam  1982). 

The  front-end  generation  program  FEGEN  allows  the  user  to  specify  a  set  of  particle  properties  such  as 
shape,  size,  contact  stiffiiesses,  dampings  and  strengths,  as  well  as  a  gradation  in  terms  of  the  relative  fraction 
of  each  property  type.  FEGEN  provides  assistance  with  choosing  the  contact  properties  based  on  the  desired 
elastic  modulus,  Poisson's  ratio  and  coelfrcient  of  restitution  of  each  particle  type.  Orientations  of  the  particles 
nity  be  specified  at  any  angle  with  respect  to  the  horizontal,  and  may  be  uniform  or  normally  distributed  with  a 
specified  standard  deviation  about  the  required  orientation.  After  the  region  bounds  and  the  desired  area  ratio 
are  specified,  the  particles  are  generated  randomly  according  to  the  gradation  and  orientation  information.  The 
particles  are  plotted  as  they  are  being  goioated.  Particles  may  be  added  or  deleted  graphically  using  the  drawing 
editcH*.  An  ASCII  data  file  is  written  which  mty  be  used  as  a  starting  input  file  to  the  simulation  program  PARSE. 

As  with  PARSE,  both  FEGEN  and  FEPARSE  were  written  and  debugged  with  BorlandC-H-  under  Microsoft 
Windows.  The  two  interface  modules  together  total  over  10,000  lines  of  C.  Because  the  interface  routines  link 
with  FEINT  which  makes  calls  to  BorlandC-specific  graphics  routines,  a  separate  translation  library  was  made 
for  Watcom  C-h-  .  When  the  interface  programs  are  compiled  with  the  WatcomC-H-  package,  the  resulting 
interface  routines  are  no  longer  confined  to  the  artifical  lower  640  KB  of  memory  limitation. 
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5.  VALIDATION 


Proper  validation  of  any  numerical  method  should  include  comparison  of  numerical  results  with  analytical 
and  physical  results.  For  particulate  systems,  few,  if  any,  analytical  solutions  are  available  where  the  system  is 
properly  treated  as  a  particulate  system.  Continuum-based  analyses  are  possible  based  on  linear  elasticity  and 
elasto-plasticity,  fix  exan^le,  but  these  would  not  offer  prq)er  validation  of  the  DEM  method  and  code  (although 
such  ccxnpariscms  would  nevertheless  be  useful  in  their  own  rigfit).  Exceptions  include  the  5 -disk  systems  subject 
to  static  equilibrium  und^  gravity  loading  (Ting  and  Corkum  1991),  or  multiple  rock  blocks  on  an  inclined  plane 
(Grabinsky  1989).  In  this  latter  case,  a  polygonal  two-dimensional  DEM  code  (ADEM)  was  used  successfully 
to  determine  the  mode  of  failure  of  a  collection  of  rock  blocks,  with  the  analytical  solution  of  the  onset  of  failure 
from  Goodman  and  Bray  (1976).  Since  analytical  solutions  for  particulate  are  rare,  the  importance  of  proper 
code  validation  with  physical  test  results  is  particularly  important, 

F(x  the  curr^t  stucty,  two  datasets  were  used  for  comparison  -  the  two  dimensional  photoelastic  rod  data  of 
Konishi,  Oda  and  Nemat-Nasser  (1983),  and  the  circular  rod  data  of  Chapuis  (1976).  Konishi  et  al  sheared 
assemblies  of  approximately  700  photoelastic  rods  of  varying  particle  angularity  (aspect  ratio)  in  order  to  study 
the  effects  of  particle  shape  <xi  mechanical  behavior.  The  detailed  physical  test  setup  and  numerical  simulations 
and  comparisons  are  described  later  in  this  section.  In  addition,  a  variety  of  particle-particle  contact  relations 
were  investigated.  These  ranged  from  contact  springs-dashpot  combinations  as  shown  in  Figure  3  where  the 
springs  are  linear  (constant  stiffhess)  to  linear  hysteretic  springs  (with  no  dashpot).  In  addition,  the  nonlinear 
ccmtact  model  proposed  by  Chang  and  Misra  (1990)  was  investigated.  However,  the  range  of  contact  forces  in 
the  current  simulaticHis  were  found  to  be  in  a  nearly  linearly  range,  and  hence  this  nonlinear  contact  relation  was 
not  explicitly  implemented  in  the  PARSE  simulation  code. 

The  next  few  sections  describe  the  linear  hysteretic  and  nonlinear  contact  laws  which  were  investigated. 


5.1  The  Hysteretic  Force-Deflection  Algorithm 

Figure  10  shows  a  linear  force-deflection  algorithm  that  includes  position  dependent  hysteresis  (Walton  and 
Braun  1986). 


Figure  10,  Linear  hysteretic  force-deflection  relationship  (Walton  and  Braun,  1986) 
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The  nonnal  force,  F„,  is  calculated  as. 


F^=K^a,loading 

(42) 

=K^(a-aQ),  unloading 

(43) 

Walton,  Brandeis  and  Cooper  (1984)  perfcraied  numaical  cqjeriments  using  an  elastic-perfectly-plastic  two 
Hiiwnginnal  finite  element  model,  NIKE2D.  The  ejqjoiments  caisisted  of  moving  a  hemisphere  into  a  rigid  wall, 
withdrawing  the  hemisphere,  and  then  moving  the  hemisphere  into  the  wall  again.  Walton  and  Braun  (1986) 
fwvtftlivt  the  resulting  force-deflection  behavior  using  a  linear  loading  curve  of  slope  K,,  and  a  linear  unloading 
curve  of  steeper  slope  Kj .  Upon  reloading,  the  normal  force  and  deformation  will  increase  linearly  with  slope 
K2  until  the  original  loading  curve  is  reached,  at  which  point,  continued  loading  will  follow  the  original  loading 
curve  with  slq)C  K,.  In  this  model,  the  area  between  the  loading  and  unloading  curves,  represents  the  ener^  lost 
to  plastic  deformation.  The  loading  and  unloading  slopes,  K]  and  Kj,  were  related  based  on  the  coefficient  of 
restitution,  e,  as  follows. 


The  programs  FEGEN  and  PARSE  were  altered  to  allow  for  hysteretic  contacts.  The  hysteretic  version  of 
FEGEN  recommends  values  fw  K,,  loading  stiffiiess,  based  on  the  relationships  previously  discussed.  For  given 
values  of  K,  and  coefficient  of  restitution,  e,  K2  is  calculated  and  stored. 


S.2  Nonlinear  Force-Deflection  Relationship 


The  nonlinear  force  deflection  relationship  of  Chang,  Misra  and  Xue  (1989)  was  also  investigated.  This 
relationship  is  a  modificaticHi  of  Hertzian  contact  theory  for  cylindrical  disks  in  contact  which  allows  for  contact 
slip.  The  normal  stiffness  k„  is; 
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In  the  preceding  equations,  G  is  the  shear  modulus  of  a  particle,  v  is  Poisson's  ratio  of  a  particle,  r  is  the 
avo’age  radius  of  the  two  particles  of  radii  ri  and  rj,  A  is  the  half-width  of  the  contact,  and  ^  is  the  normal  contact 
force.  Fa- two  ellipses  in  contact  (instead  of  disks),  modifications  were  made  to  these  equations.  The  average 
radius  r  used  in  the  solution  for  ellipses  in  contact  was  obtained  by  first  averaging  the  major  and  minor  axes  of 
each  ellipse,  and  then  taking  the  average  for  the  two  ellipses. 


In  this  nonlinear  contact  model,  the  normal  force  increases  nonlinearly  as  the  deformation  at  the  contact 
Similarly,  the  normal  contact  stiffness  varies,  and  is  defined  as  the  instantaneous  slope  for  any  two 
values  of  normal  force,  deformation,  delta,  as  illustrated  in  Figure  1 1. 
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which  simplifies  to; 

A  J-)  fAdA --  fAlnAdA 

r  2r  J  rJ 

Integrating  this  equation  yields; 

A=(i2it-±)i^  \lnA)-^A 

r  2r  2  2r  2 


(52) 


(53) 

(54) 


Using  equaticn  (54)  to  calculate  the  half-width.  A,  and  the  deformation,  A,  for  input  values  of  normal  force, 
Fn,  widi  appropriate  values  of  Young’s  modulus  and  Poisson’s  ratio  for  rubber  results  in  the  nonlinear  contact 
relation  plotted  in  Figure  12  below. 


Figure  12.  Nonlinear  force-deflection  relationship  from  Chang  et  al  (1989). 

To  implement  the  nonlinear  force-deflection  relationship,  it  was  necessary  to  obtain  a  simple  equation  that 
would  relate  the  particle  overlap  at  a  contact  to  the  normal  force.  A  least  squares  curve  fit  of  the  data  generated 
in  Figure  12  was  carried  out  using  a  power  flinction.  This  resulted  in  the  following  equation  fw  normal  force  as 
a  fimction  of  normal  contact  overlap  delta; 

F,  =  aj  A  +  ajA^  +  ajA^  (55) 

where; 

a,  =4786.786;  a2=  12330798;  a3=-6.8e9 
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Since  the  overlap  at  contacts  in  the  DEM  simulations  is  small  compared  to  the  radii  of  the  particles,  a  linear 
regressi(xi  analysis  of  the  pcMticm  of  the  curve  below  a  normalized  deformation  of  0.02  was  deemed  appropriate. 
The  purpose  of  this  analysis  was  to  determine  how  close  the  data  was  to  being  linear  in  this  region  of  interest. 
The  result  of  the  linear  regressm  analysis  are  plotted  together  with  the  modified  Chang  force  deflection  relation 
in  Figure  13.  The  linear  regression  analysis  of  these  values  yielded  a  slope,  k^,  of  5783  kN/m,  and  an  value 
of  0.985.  In  this  small  contact  deformation  range,  the  normal  force-deflection  relation  is  nearly  linear. 


Figure  13.  Best  least  squares  linear  force  contact  relation  compared  with  modified  Chang  relation. 


5.3  Chapuis  (1976)  Test  Data 


Chqjuis  (1976)  perfOTined  biaxial  compression  tests  on  assemblies  of  400  to  600  rubber  cylinders,  each  10 
mm  in  length.  The  cylinders  were  of  four  differoit  sizes,  with  diameters  of  1 9.05  mm,  22.38  mm,  25.65  mm,  and 
38.76  mm,  with  gradation  proportions  of  4;3:2: 1  respectively.  The  two  largest  cylinders  were  nearly  circular 
in  cross  section,  while  the  smaller  diameter  particles  were  reported  to  be  between  2  to  5  percent  out  of  round. 
The  resulting  aspect  ratio  of  the  smaller  particles  was  approximately  1.03. 

Chapuis  performed  several  tests,  including  a  simple  shear  test  of  a  dense  assembly,  biaxial  tests  on  both  a 
dense  and  a  loose  assembly,  and  a  biaxial  unloading  test.  As  part  of  DEM  model  validation,  only  the  dense  and 
loose  biaxial  tests  and  the  unloading  test  were  attonpted.  These  results  were  reported  in  detail  in  Rowell  ( 1 994), 
and  are  summarized  in  this  Report. 

From  the  reported  data,  the  initial  isotropic  confining  stress  prior  to  shear  was  approximately  3.6  Ib/in  cm, 
or  63  kN/m^.  The  biaxial  compressicxi  testing  apparatus  used  to  shear  the  assemblies  is  shown  in  Figure  14.  The 
size  of  the  ^paratus  in  Figure  14  was  measur^  to  arrive  at  the  initial  height  to  width  ratio  of  the  sample.  This 
infixmaticai  was  then  combined  with  gradaticxi  and  void  ratio  infcxmation  to  arrive  at  an  approximate  total  sample 
area  of 0.3267  m^  and  an  initial  sample  size  of  0.484  m  width  by  0.675  m  height. 

The  fi'ont-end  graphical  interface  program,  FEGEN,  was  used  to  create  dense  and  loose  samples  with  the 
same  gradation  and  niunber  of  particles  as  the  Chapuis  tests.  Figure  15  represents  one  sample  that  was 
compressed  to  form  a  dense  assembly  of  particles. 

Table  2  contains  the  particle  size  infcxmation,  gradaticxi  information,  and  mcxlel  input  parameters  for  the  first 
attempt  to  simulate  the  Chapuis  tests.  One  of  the  difficulties  encountered  in  modeling  these  tests  was  that  the 
type  and  properties  of  the  material  being  tested  were  not  reported.  (  It  has  since  been  determined  by  direct 
ccxnmunicadcxi  that  the  rods  woe  made  from  hard  rubber).  Since  the  confining  stress  applied  for  these  tests  was 
sevoal  orders  of  magnitude  greater  than  that  applied  for  the  Konishi  et  al  tests,  this  suggested  that  a  much  stiffer 
material  was  used  in  the  Chapuis  tests.  The  stifeess  values  used  in  Table  2  resulted  in  large  particle  overlapping 
and  an  overall  material  response  that  was  much  tcx)  soft.  Subsequent  simulations  included  (higher)  stiffiiess 
values  which  wcxild  more  accurately  model  the  physical  test  results.  For  example,  the  interface  program  FEGEN 
recommended  a  value  of  k„=  70,000,000  kN/m  for  disks  made  from  granite.  Simulations  using  higher  normal 
contact  produced  excellent  agreement  with  the  physical  results.  The  stiffness  value  used  in  these  successful 
simulations  was  k„=  2,400,000  kN/m.  Also,  the  physical  interparticle  sliding  fnction  reported  for  the  ovals  was 
16°,  not  the  10°  u^  in  the  simulation.  According  to  Chapuis,  the  reported  1 0°  fiiction  angle  includes  particle 
rolling,  while  the  reported  16°  fnction  angle  was  for  pure  sliding. 

Figure  16  presents  the  results  from  attempts  to  model  the  Chapuis  tests.  Two  different  samples  were  created 
and  tested.  Both  samples  were  tested  with  the  same  model  parameters.  These  simulations  were  conducted  with 
k„=  2,400,000  kN/m.  The  maximum  friction  angle  obtained  during  the  physical  testing  of  the  dense  sample  was 
25°.  The  DEM  simulations  yielded  a  maximum  fnction  angle  of  24  ,  which  is  in  excellent  agreement.  The 
volume  change  behavior  of  the  DEM  simulated  test  also  compared  favorably  with  the  Chapuis  tests.  This  data 
is  presented  in  Figure  17. 
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Figure  14.  Testing  apparatus  used  in  Chapuis  (1976)  tests. 


Figure  15.  DEM  sample  for  dense  Chapuis  test. 
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K^tanca  befOK  contact  list  brcdEca 
accumulatof  toicmtces  for  deflection  and  rotatiM 
coefBcieiit  of  restitution 


TABLE  2.  Simulation  parameters  for  an  initial  simulation  of  Chapuis’  tests. 
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Figure  16.  Validation  results  for  dense  Chapuis  (1976)  disk  tests:  normalized  deviator  stress  vs. 
axial  strain. 
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Figure  17.  DEM  validation  results  of  for  dense  Chapuis  test;  volumetric  strain  vs.  axial  strain. 

Efifoits  were  also  made  to  simulate  a  loose  biaxial  shear  test,  and  an  unloading  test  performed  by  Chapuis. 
The  attonpts  to  simulate  the  loose  test  were  unsuccessful,  due  to  difficulties  associated  with  obtaining  a  sample 
at  the  same  void  ratio  and  initial  stress  as  the  Chapuis  loose  sample.  The  same  problems  were  encoimtered  in 
modeling  the  Ch^uis  unloading  test.  However  an  i^oading  test  was  perfixmed  at  an  initial  void  ratio  different 
than  the  value  rqxirted  by  Ch^uis,  and  the  results  are  shown  in  Figure  18  for  a  qualitative  compariscm.  As  this 
Figure  suggests,  the  behavior  of  the  sample  is  qualitatively  comparable  to  the  unloading  test  performed  by 
Chapuis,  with  the  sample  rebounding  to  a  denser  state.  It  is  believed  that  better  agreement  is  possible  if  a 
numerical  sample  could  be  created  with  the  same  initial  density  as  the  Chapuis  unloading  test. 

Overall,  the  DEM  results  woe  similar,  both  qualitatively  and  quantitatively,  to  the  physical  tests  performed 
by  Chq)uis  ai  elliptical  particles.  The  measured  peak  friction  angle,  4)n^,  from  the  DEM  simulation  of  a  dense 
sample  was  24%  ccmipared  to  25°  achieved  during  the  physical  tests.  Volumetric  strain  data  also  compared 
favorably.  Although  the  Chapuis  tests  were  conducted  without  complete  information  available  regarding  the 
material  being  tested,  it  is  conjectured  that  modeling  of  the  Chapuis  tests  was  successful  because  the  particle 
shape  and  gradation  matched  in  the  physical  and  numerical  tests.  Additionally,  the  stiffiress  value  used  in  the 
simulation  was  large  enough  to  guarantee  small  particle  deformations,  which  is  the  key  to  achieving  proper 
mechanical  behavior  in  these  numerical  simulations. 
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Figure  18.  Results  of  attempted  simulation  of  Chapuis(1976)  unloading  test. 
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5.4  Konishi,  Oda  and  Nemat-Nasser  (1983)  Tests 

Koiishi  etal.  (1983)  pofcxmed  physical  tests  on  assemblages  of  approximately  700  photoelastic  oval  rods 
of  vaiying  aspect  ratio,  subjected  to  biaxial  shear.  This  chapter  describes  the  physical  test  setup  and  results  of 
these  tests.  Etetailed  photographic  prints  of  one  of  the  physical  tests  were  obtained  and  digitized  for  use  in  one 
of  die  simulations.  Model  input  parametos  were  varied  extensively  to  assess  the  effect  on  the  simulation  results. 

The  oval  {diotoelastic  rods  had  common  length  of  19  mm.  Two  different  cross  sectional  shapes  were  used; 
“Oval  I”  and  “Oval  IF,  with  1 . 1 : 1  and  1 .4;  1  particle  aspect  ratios.  Samples  were  formed  from  three  different 
sized  particles  for  eadi  oval  aoss  sectional  shape.  Table  3  summarizes  the  particle  dimensions  used.  Two  sets 
of  experiments  were  performed.  In  one  set,  the  particles  were  lubricated  with  talcum  powder  that  reportedly 
prodiKed  an  estimated  interparticle  friction  angle  of  26  ° .  The  other  set  consisted  of  unlubricated  particles  with 
an  estimated  interparticle  friction  angle  of  52^“.  The  rods  were  cast  from  polyurethane  rubber  with  a  reported 
100%  modulus  of  80  kgf/an^.  Particle  cdiesioi  was  reportedly  zm)  for  all  tests.  Figure  19  illustrates  the  biaxial 
testing  rqrparatus  used  to  test  the  samples.  Initial  dimensions  of  the  loading  frame  were  330  rtun  height  and  207 
mm  wid&  The  particles  were  stadced  by  hand  in  the  loading  frame,  which  could  be  tilted  to  create  samples  with 
varying  bedding  angle.  Samples  with  bedding  angles  of  0°,  30°,  60°,  and  90°  were  tested. 


Designation 

a^ 

Diameters  a^  (mm) 

Oval  I 

1. 1 

Large  14.8/13.4 

Medium  9.9/8.9 

Small  6.3/5.7 

Oval  II 

1.4 

Large  16.0/11.3 

Medium  10.7/7.4 

Small  7. 1/4.9 

Table  3.  Description  of  particle  sizes  and  shapes  used  by  Konishi  et  al  (1983). 

A  constant  face  of  0.45  kgf  (0.70  kPa),  was  applied  laterally  to  the  samples,  and  vertical  compression  was 
appUed  incrementally  by  lowering  the  upper  part  of  the  loading  frame.  Displacements  were  monitored  by  dial 
gauges,  and  load  transducers  monitored  applied  loads.  The  side  platens  of  the  loading  frame  were  coated  with 
Teflon  to  reduce  wall  friction,  and  the  applied  horizcmtal  and  vertical  stresses  were  assumed  to  be  principal 
stresses.  Throughout  the  loading,  photoelastic  photographs  were  taken  to  record  force  chain  formation.  Table 
4  presents  the  results  from  the  tests  conducted  by  Konishi  et  al.  (1983). 

Samples  that  matched  the  general  gradation,  aspect  ratio,  number  of  particles,  and  bedding  angles  of 
Konishi’s  tests  were  generated  using  the  front-end  graphical  interface  program  FEGEN.  The  results  of  these 
initial  validation  tests  have  been  reported  by  Ting  et  al  (1993),  and  exhibited  qualitatively  similar  behavior  to 
Konishi's  tests.  However,  the  peak  stress  ratios  reported  by  Konishi  et  al  were  4.9  for  Oval  I  samples  of  26° 
intonal  frictiem,  0°  bedding  angle,  and  1 1 .8  for  Oval  II  samples  of  26°  internal  friction,  0°  bedding  angle.  These 
values  of  peak  stress  ratio,  and  consequoitly  peak  friction  angle,  were  coisistently  higher  than  the  resulting  values 
of  DEM  simulations.  Figure  20  presents  the  results  of  the  initial  DEM  validation  tests,  plotted  as  peak  friction 
angle  versus  aspect  ratio.  Additional  tests  are  presented  later  in  this  Report,  as  well  as  Meachum  (1994). 
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peak  friciion  angle  0  (®) 


1  loading  fraa* 

2  load  tranaducar 

3  no  tor  for  vertical  loading 

4  dial  gauge 

5  LVDT'a 

6  weights  for  lateral  loading 

7  tilting  axes 


Figure  19.  Schematic  of  biaxial  testing  apparatus 
from  Konishi,  Oda  and  Nemat-Nasser  (1983). 


Figure  20.  Peak  angle  of  internal  friction  measured  in  initial  DEM  simulations  of  biaxial 
shear  tests  as  a  function  of  aspect  ratio  (Ting  et  al  1993). 


Cross  Friction  Angle  Bedding  Angle  Void  Ratio  e  Peak  Stress  Internal 
Section  <!>„  0  Ratio  (ol/o2)  Friction 

Angle  <|)^ 


0“  0.177  14.4  60° 

30°  0.154  11.4  57° 


OVAL  II  52° - 

60°  0.159  8.3  52° 


OVAL  II  26° - 

60°  0.155  4.3  39° 


_ I _ I  90°  I  —  I  5.5  I  44° 

Table  4.  Results  of  physical  tests  on  oval  rods  performed  by  Konishi  et  al  (1983). 


For  the  sanples  described  above,  Koiishi  et  al  rqxxt  peak  friction  angles  of  42°  for  Oval  I,  and  58°  for  Oval 
II  for  the  0°  bedding  angle  samples.  The  initial  sequence  of  DEM  simulations  reported  by  Ting  et  al  (1993), 
yielded  peak  fricticxi  angle  values  of  3 1°  fw  Oval  I,  and  37°  for  Oval  II.  It  was  therefore  necessary  to  investigate 
the  reasons  for  the  discrepancy  in  the  macroscopic  strength  behavior  of  the  DEM  samples. 

To  bettCT  model  the  system  of  particles  tested  by  Konishi  et  al,  photoelastic  photos  of  the  Oval  II  sample 
with  0°  bedding  angle,  were  obtained  from  the  authors  and  digitized.  The  Oval  II  sample  that  was  obtained  and 
digitized  is  dqricted  in  Figure  21.  The  digitized  data  was  used  to  create  a  numerical  sample  that  would  accurately 
represoit  the  exact  particle  locations  and  orientations  of  Konishi's  Oval  II  sample.  The  sample  was  created  with 
the  program  FEGEN,  and  is  represented  in  Figure  22. 
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Figure  21 .  Photo  of  Konishi  et  a/’s  Oval  II  sample  prior  to  shear. 


Figure  22.  Digitized  version  of  Konishi  et  al 's  Oval  II  sample. 
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Figure  23.  Comparison  of  initial  DEM  simulation  of  biaxial  shear  using  digitized  Oval  II 
sample. 
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5.5  Variation  of  Contact  Properties 

Althou^  careful  coosido'aticMi  was  given  to  the  pr<^)erties  and  simulation  parameters  used  in  the  initial  DEM 
tests,  additional  simulations  with  varying  contact  properties  were  conducted  on  the  digitized  Oval  II  sample. 
These  tests  were  conducted  to  assess  the  effect  of  variation  of  simulation  parameters,  and  to  determine  if  the 
DEM  simulaticxis  would  more  closely  model  the  Konishi  tests.  Extensive  testing  was  performed  as  part  of  this 
validation  procedure,  consisting  of  contact  stiffiiess  variation,  the  addition  of  contact  cohesion,  variation  of 
contact  friction,  variation  of  contact  damping,  the  addition  of  wall  friction,  testing  with  gravitational  effects 
implemented,  and  simulations  with  various  combinations  of  these  parameters  invoked.  The  effect  of  changing 
the  timestep  of  integraticm  was  studied,  as  well  as  the  effect  of  shearing  the  samples  at  different  rates.  These  test 
results  are  describ^  in  full  in  Rowell  (1994),  and  only  the  major  findings  are  reported  here. 

5.5.1  Variation  of  Linear  Contact  Spring  Stiffness 

To  determine  the  effect  of  contact  stiffiiess  variation  on  the  DEM  simulation  results,  the  linear  contact 
stiffness  k„  was  varied  from  300  kN/m^  to  10,000  kN/irf  while  all  model  parameters  not  affected  by  contact 
stiffness  remained  constant.  Candle  and  Strack  showed  that  when  the  ratio  of  tangential  stiffness  to  normal 
stiffiiess,  k^  was  varied  frcMn  2/3  to  1,  there  was  little  change  in  mechanical  behavior.  Of  the  nine  simulations 
that  were  perfcxmed  as  part  of  this  testing  sequence,  eigjit  simulations  used  k,  equal  to  2/3  k^.  The  remaining  test 
was  run  with  k^  equal  to  k,.  Figure  24  illustrates  the  effect  of  contact  stiffiiess  variation  on  the  DEM  results. 


Figure  24.  Effect  of  varying  linear  normal  contact  stiffiiess  k„  on  DEM  simulation  results. 


37 


5.5.2  Alternate  Contact  Stiffness  Model 


DEM  simulations  were  also  conducted  using  the  linear  hysteretic  contact  model  proposed  by  Walton  and 
Braun  (1986).  For  a  specified  coefficient  of  restitution  of  0.5  and  a  loading  stiffiiess  of  K,=  4500  kN/m^,  the 
unloading  stiffiiess,  Kj  was  calculated  as  18,000  kN/m^.  All  other  simulation  parameters  and  contact  properties 
remained  unchanged  from  the  linear  contact  parameters  except  for  the  contact  damping,  which  is  zero  in  the 
hysteretic  model.  Figure  25  plots  die  nomalized  deviator  stress  versus  axial  strain  for  the  linear  spring-dashpot 
and  linear  hysteretic  spring  tests,  and  shows  that  these  two  contact  models  exhibit  similar  behavior. 


vertical  strain 

Figure  25.  Comparison  of  linear  spring-dashpot  and  linear  hysteretic  spring  ccmtact 
models  for  biaxial  shear  of  digitized  Oval  II  sample. 


5.5.3  Variation  of  Contact  Cohesion 

Although  Kcaiishi  et  al  repcxted  that  the  photoelastic  material  used  in  the  physical  tests  possessed  no  contact 
cdiesion,  it  is  possible  that  a  small  amoimt  of  cohesion  existed  at  the  particle  contacts.  To  assess  the  effect  of 
this  cdiesicxi,  DEM  simulaticms  were  carried  out  where  a  small  cohesion  was  introduced.  The  shear  force,  T,  at 
a  contact  is  related  to  the  normal  force,  N,  at  a  contact  by  the  relationship, 

T  =  „+Ntan<|>„ 

where  c..  is  the  contact  cohesion,  and  <|\,  is  the  interparticle  friction  angle.  This  relationship  was  used  with  the 
results  firxn  previous  DEM  tests  to  arrive  at  realistic  cohesicm  values  for  input  into  the  model.  The  DEM  results 
from  this  series  of  tests  are  represented  in  Figure  26,  and  showed  no  appreciable  difference  over  the  range  of 
cohesion  values  tested. 
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Figure  26.  Effect  of  varying  contact  cohesion  on  DEM  simulation  results  for  digitized  Oval 
II  sample.. 


5.5.4  Variation  of  Interparticle  Friction  Angle,  <|>, 

The  average  interparticle  friction  angle  reported  by  Konishi  et  al.  (1983)  for  their  particles  lubricated  with 
talcum  powder  is  26°.  How  this  friction  angle  was  measured  was  not  described  in  the  available  literature.  In 
addition,  one  can  envisage  how  this  friction  angle  value  could  change  during  testing  as  a  result  of  particles  sliding 
past  each  other,  resulting  in  an  uneven  distributiai  of  talcum  powder.  Since  this  friction  angle  value  is  an  average 
value,  and  because  thoe  is  s(xne  uncertainty  related  to  this  value,  DEM  simulations  were  performed  on  samples 
with  varying  inteiparticle  friction  angle  while  all  other  parameters  remained  constant.  The  results  of  these  tests 
are  plotted  in  Figure  27. 


5.5.5  Effect  of  Contact  Damping,  c. 

The  normal  contact  damping  value,  Cn ,  can  be  related  explicitly  to  the  contact  stiffness,  masses  of  the 
particles  in  contact,  and  the  coefficient  of  restitution  of  the  particles,  according  to  equation  [35]  (  Corkum  and 
Ting  1986).  SiiKe  the  value  of  the  coefficient  of  restitution  of  the  photoelastic  material  used  in  the  physical  tests 
is  not  reported  in  the  available  literature,  a  value  had  to  be  selected  arbitrarily.  For  rubber,  a  realistic  value  for 
the  coefiBcient  of  restitutiixi  was  selected  as  e  =  0.5,  based  on  typical  known  values.  All  of  the  simulations  in  this 
study  used  values  ofe  equal  to  0.5.  To  learn  if  higher  damping  would  significantly  influence  simulation  results, 
a  test  with  e  =  0.25  was  conducted.  The  resulting  data  was  plotted  in  Figure  28. 
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5.6  Discussion  of  Validation  Testing  and  Konishieta/  Data 

The  value  of  47°  obtained  in  the  DEM  simulation  is  still  much  lower  than  the  reported  value  of  58°  obtained 
in  the  physical  tests.  Konishi  reported  a  maximum  normalized  deviator  stress  of  1 1.8,  which  was  achieved  at 
^proximately  1.5%  axial  strain,  fex*  the  Oval  n  sample  with  an  interparticle  friction  of  26°  and  0°  bedding  angle. 
The  DEM  results  foe  this  test  yielded  a  maximum  normalized  deviator  stress  of  approximately  5.5  at  4%  strain. 

Sev^inl  possibilities  exist  that  would  e?q)lain  the  differences  between  the  numerical  and  physical  test  results. 
One  possibility  is  that  the  linear  ccxitacts  used  in  the  DEM  model  does  not  adequately  reflect  the  physical  contact 
behavior.  This  possibility  has  been  investigated  through  the  study  of  a  hysteretic  contact  relationship,  and 
through  the  analysis  of  a  nonlinear  contact  relationship.  Comparisons  between  simulation  results  incorporating 
Waltcxi's  hysteretic  contact  model  and  the  simple  linear  contact  model  currently  xised  in  the  simulation  program 
PARSE,  reveal  litde  difference  in  the  observed  macroscopic  mechanical  behavior  of  the  samples  tested.  Figure 
25  suggests  that  both  the  linear  hysteretic  spring  and  the  linear  spring-dashpot  contact  models  yield  results 
coisistoit  with  each  other  and  with  typical  soil  behavior.  This  research  finding  suggests  that  although  these  two 
contact  models  are  modeling  the  force-deflection  behavior  and  energy  dissipation  at  particle  to  particle  and 
particle  to  wall  contacts  differently,  the  overall  mechanical  response  yielded  by  both  contact  models  is 
comparable. 

Analysis  of  the  ncxilinear  Hertzian  contact  model  drived  by  Chang  and  Misra  (1990)  suggests  that  for  small 
ccxitact  defcxmaticMis,  the  relationship  between  normal  contact  force  and  normal  deformation  very  nearly  linear 
Since  it  is  necessary  to  use  stiffiiess  values  in  the  DEM  model  that  will  result  in  small  deformation  or  overlap  of 
particles  relative  to  particle  size,  analysis  of  the  nonlinear  contact  model  focused  on  contact  deformations  in  the 
range  of  iqj  to  2%  of  the  avCTage  particle  radius.  It  was  found  that  in  this  range  of  small  contact  deformation,  the 
resulting  force-deformation  relationship  is  nearly  linear,  with  a  resultant  value  of  0.985.  Although  Hertzian 
contact  theory  suggests  that  the  force-deformation  relationship  of  solids  in  contact  is  nonlinear  in  nature,  it  is 
unnecessary  to  model  this  relationship  exactly  in  the  DEM  model  for  the  current  low  stress  level  simulations. 

Based  on  these  findings,  and  the  findings  concerning  comparisons  with  the  hysteretic  contact  model,  it 
appears  that  the  DEM  model  with  the  simple  linear  caitact  relationship  is  sufficient  to  properly  model  the  overall 
mechanical  behavior  of  particulate  systems.  The  key  to  the  proper  selection  of  normal  contact  stiffiiess  for  the 
simple  linear  DEM  model  is  that  deformation  at  particle  contacts  must  be  kept  small  relative  to  the  size  of  the 
particle.  It  is  unlikely  that  the  force-deflection  relationship  is  responsible  for  the  differences  between  the 
niunerical  and  physical  test  results. 

The  extensive  validaticxi  effext  ccxisisting  of  variation  of  DEM  model  input  parameters  also  failed  to  provide 
a  better  e^proximaticxi  of  Kcxiishi  et  a/’s  physical  test  results.  Only  small  changes  in  the  macroscopic  behavior 
were  observed  through  this  effort.  Therefore,  it  is  highly  unlikely  that  the  reason  for  the  discrepancy  between 
physical  and  numerical  test  results  lies  in  the  improper  selection  of  model  input  parameters. 

It  is  conjectured  that  the  reason  for  the  discrepancy  in  test  results  lies  in  the  fact  that  the  photoelastic  rods 
sheared  in  the  Konishi  et  al  tests  arc  oval  and  not  elliptical  in  shape.  Figure  29  presents  an  enlarged  view  of  the 
lower  left  comer  of  the  Oval  II  photoelastic  photo  taken  at  stage  1,  provided  by  Prof.  Oda.  A  close  inspection 
of  the  particle  shape  reveals  that  the  Oval  II  rods  have  a  sigmficantly  large  flat  surface  parallel  to  the  major  axis. 
In  addition,  for  the  0°  bedding  sanq)le,  these  flat  surfaces  arc  aligned  horizontally,  which  creates  a  situation  where 
many  particles  are  in  contact  over  these  large  flat  regions.  When  the  sample  is  loaded  vertically,  the  load  can  be 
transferred  directly  through  chains  of  these  stacked  flat  rods,  which  are  very  stable. 
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Figure  29.  Enlarged  view  of  lower  left  comer  of 
Konishi  et  al  (1983)  Oval  11  sample  at  stage  1. 

As  further  siqipcxtof  this  theoiy,  simulations  at  bedding  angles  greater  than  0°,  specifically  the  60°  bedding 
samples,  yielded  Action  angles  much  closer  to  the  Konishi  test  results.  Konishi  reports  a  friction  angle, 
of  39°,  fOT  the  Oval  H,  26  interparticle  friction,  60  bedding  test  sequence.  A  DEM  simulation  of  this  test 
sequence  yielded  a  fiicticxi  angle,  of  37°(Meachum  1994)  For  these  samples,  the  effect  of  the  flat  surfaces 
is  less  praiounced,  as  diese  flats  are  not  aligned  with  the  loading  platen,  resulting  in  a  less  stable  sample  fabric. 

To  quantify  the  particle  sh^  differences  that  exist  between  the  photoelastic  oval  rods  and  the  ellipses  used 
to  nnodel  the  ovals,  it  is  useful  to  look  at  the  contact  normal  distributions  of  each  system  of  particles.  Figure  30 
plots  the  distributiai  of  contact  normals  for  the  Oval  II,  26°  interparticle  friction,  0°  bedding  angle  Konishi  et  al 
test  at  the  initial  coiditicxi  pricH*  to  shear.  Figure  3 1  plots  the  contact  normal  distribution  at  the  peak  stress  state. 

Fa*  con:q)aris(xi,  the  contact  ncxmal  distributicxis  for  the  initial  and  peak  stress  states  of  the  DEM  simulation 
are  presented  in  Figures  32  and  33.  As  expected,  the  contact  normals  in  the  initial  state  are  generally  directed 
perpendicular  to  the  bedding  plane  (ie  vertically).  However,  for  the  DEM  simulation,  the  contact  normals  are 
more  widely  distributed  about  Ae  90°  direction.  Since  the  particle  locations  and  orientations  are  identical  at  the 
initial  condition  fa*  both  the  numerical  and  physical  tests,  ^e  difference  in  the  contact  normal  distribution  must 
be  due  to  the  particle  shape  differences. 

To  help  further  explain  the  effect  of  the  flatter  oval  particles  in  the  physical  tests,  the  physical  and  DEM 
samplfis  are  shown  at  pe^  stress  state  in  Figures  34  and  35.  Inspection  of  these  Figures  shows  that  the  contact 
tKxmals  are  aligned  nearly  vertically  in  the  physical  tests.  In  addition,  it  is  clear  that  a  very  stable  soil  fabric  has 
bear  formed  with  the  stacked  oval  rods.  Figure  35  shows  that  the  DEM  ellipse-based  sample  possesses  contact 
normals  which  are  more  widely  distributed  and  in  all  likelihood  a  soil  fabric  that  is  less  stable  than  the  system 
of  ovals. 


Figure  30.  Distribution  of  contact  normals  from  Konishi  et  al  (1983)  Oval  II  test 
prior  to  shear. 


** 


Figure  31.  Distribution  of  contact  normals  for  Konishi  et  al  Oval  II  test  at  peak 
stress  state. 

Further  compariscm  with  plQ^ical,  two  dimoisioial  test  data  are  required  to  fully  validate  the  numerical  code. 
Additional  data  are  presently  being  collected  on  the  interfacial  shear  of  photoelastic  rods  and  a  variety  of  smooth 
and  rough  surfaces  (Paikowsky  and  Xi  1995).  As  part  of  this  research  program,  contact  normal  and  tangential 
forces  are  being  computed  based  on  photoelastic  fringe  patterns.  Rods  which  are  carefully  machined  to  be 
circular,  and  also  ellipse-shaped  are  being  tested  in  these  shear  tests.  These  physical  disk  and  eventually  ellipse- 
based  physical  test  results  will  be  used  to  help  validate  the  PARSE  numerical  code  (Paikowsky,  Ting,  Xi  and 
Mischel,  in  press  1996). 
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Figure  32.  Contact  normal  distribution  of  1.4;  1  aspect  ratio  DEM  sample  prior  to  shear. 


Figure  33.  Contact  normals  distribution  of  1.4;  1  aspect  ratio  DEM  Simulation  at  peak  stress. 
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Figure  34.  Photo  of  Konishi  et  al.  ’s  Oval  II  horizontally  bedded 
sample  at  peak  stress  state. 
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Figure  35.  Horizontally  bedded  1.4:1  : 1  aspect  ratio  sample  at  peak  stress  state  based  on  DEM 
simulation  of  initial  digitized  Oval  II  sample.. 


5.7  Comparison  with  Other  DEM  Results 

5.7.1  Numerical  Results  from  Rothenburg  and  Bathurst  (1991) 

Several  researchers  have  used  ellipse-based  DEM  models  to  stucfy  particle  shape  effects  on  the  mechamcal 
behavia-  of  particulate  systems.  Rothenburg  and  Bathurst  (1991, 1992)  performed  biaxial  compression  tests  on 
assemblies  of  1000  elliptical  shaped  particles  of  varying  density  and  eccentricity.  Samples  were  formed  by 
growing  the  particles  from  the  same  initial  assembly.  This  was  accomplished  by  selecting  a  target  eccentricity, 
which  was  inoeased  incrementally  while  maintaining  hydrostatic  boundary  conditions.  The  assemblies  created 
in  this  manner  were  dense,  with  no  preferred  bedding  angle.  Since  no  particle  size  is  specified  in  this  generation 
grhpmf!^  particle  gjypis  vary.  A  simple  linear  contact  model,  similar  to  that  used  by  Cundall  and  Strack  and  in  this 
study,  was  used  to  govern  contact  behavior.  For  all  tests  conducted,  the  interparticle  fiiction  angle  was  26.6°. 

Rothenburg  and  Bathurst  express  the  angularity  of  the  ellipses  in  their  stu^  as  particle  eccentricity,  e,  which 
is  defined  as, 

e= — -  (56) 

a+b 

where  a  and  b  are  the  values  of  the  majcM-  and  minor  axis  lengths  respectively.  A  plot  of  the  peak  friction  angle 
vs.  eccentricity  obtained  in  the  tests  performed  by  Rothenburg  and  Bathurst  (1991)  is  presented  in  Figure  36. 

Figure  36  shows  that  the  maximum  fiiction  angle  obtained  by  Rothenburg  and  Bathurst  was  42°  at  an 
eccentricity  of  approximately  0.15,  which  is  equivalent  to  an  aspect  ratio  of  1.35.  As  particle  eccentricity 
increases  beyond  a  value  of  0.15,  the  shear  strength  decreases. 
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Figure  36.  Peak  friction  angle  vs.  eccentricity  from  Rothenburg  and 
Bathurst  (1991) 
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5.7.2  Numerical  Results  from  Ng  (1992) 


Ng  (1-992)  performed  numerical  simulations  on  assemblies  of  elliptical  particles,  using  a  DEM  program 
ELLIPSE2.  As  with  the  current  study,  simulaticxis  were  performed  on  assemblies  of  varying  particle  aspect  ratio 
to  assess  the  effect  of  particle  angularity  on  the  mechanical  behavior  of  the  specimens. 

The  simulations  were  conducted  using  input  parameters  for  quartz  particles  with  a  shear  modulus  G  of  2.9 
X  10^  MPa,  Poissoi's  ratio  0.  IS,  and  a  particle  to  particle  sliding  friction  angle  of  26.6°.  The  specimens  consisted 
of  particles  of  (xie  size  only  for  each  sample.  Samples  were  generated  with  particle  aspect  ratios  of  1.00,  1.25, 
l.Sl  and  1.67.  Sample  generation  was  accomplished  with  a  depositional  particle  generation  scheme,  which 
randomly  generated  samples  by  introducing  a  gravity  field  while  particles  were  being  rained  from  the  top  of  the 
sanq>le.  The  q)ecimens  were  consolidated  isotropically,  following  sample  generation,  in  increments  of  40  kPa, 
iq)  to  iqjproximately  160  kPa  total  consolidation  pressure.  Prior  to  the  shearing  of  the  samples,  Ng  reports  void 
ratios  varying  between  0.254  for  the  sample  with  circular  elements,  to  approximately  0.23  for  the  samples  of 
aspect  ratio  1.51  and  1.67. 

Biaxial  coiiq)ressicm  tests  ware  ccmducted  on  three  of  the  samples.  The  sample  with  particles  of  aspect  ratio 
1 .67  was  sheared  to  12%  axial  strain.  The  samples  of  aspect  ratio  1 .00  and  1 .25  were  sheared  to  only  2.5%  axial 
strain,  due  to  computational  time  constraints.  Ng  fcxuid  that  the  particles  with  larger  aspect  ratios  produced  larger 
friction  angles,  a  result  that  is  consistent  with  the  results  produced  by  the  PARSE  program.  The  friction  angles 
yielded  by  the  ELL1PSE2  program  were  low,  however,  compared  to  typical  friction  angles  for  medium  dense 
natural  sands.  Figure  37  illustrates  the  results  obtained  by  Ng  ( 1 992)  for  the  three  biaxial  tests  performed.  Note 
that  in  Figures  37  and  38,  the  deviator  stress  is  quantified  using  the  mobilized  friction  angle,  which  is  just: 

mobilized  <|)  =  sin‘*(^)  (56) 

where  p  and  q  are  the  mean  and  maximum  deviator  stress,  respectively. 

Figure  39  plots  the  simulation  results  on  ellipse  assemblages  of  Rothenburg  and  Bathurst  (1991),  Ng  1991) 
and  DEM  program  PARSE.  The  fact  that  the  assemblies  created  by  Rothenburg  and  Bathurst  were  dense,  with 
no  preferred  bedding  angle  accounts  fix  the  increased  strength  evident  in  this  comparison.  The  assemblies  tested 
with  the  PARSE  program  were  horizontally  bedded,  medium  dorse,  and  consisted  of  three  different  particle  sizes. 
The  comparison  of  data  is  quite  good  in  light  of  these  differences. 

The  assemblies  created  by  Ng  were  loose  ccxiqrared  to  the  sanqrles  tested  with  the  PARSE  program  and  those 
tested  by  Rothenburg  and  Ba^ursL  As  expected,  the  resulting  shear  strength  of  the  Ng  samples  is  less  than  that 
yielded  by  the  other  two  programs.  Additionally,  Ng  tested  samples  consisting  of  one  particle  size  as  opposed 
to  the  graded  samples  of  the  other  two  researchers.  In  general,  the  PARSE  results  compare  favcxably  with  the 
other  DEM  model  results. 
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Fieure  38.  Mobilized  friction  angle  vs.  axial  strain  for  1 .67: 1  aspect  ratio  sample  from  Ng 
(1991). 
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Figure  39.  Plot  of  maximum  angle  of  internal  friction  vs.  particle  eccentricity  e  for 
three  DEM  models. 
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6.  SAMPLE  FORMATION  AND  TEST  SIMULATIONS 


Tlie  next  few  Sections  deal  with  the  use  of  the  ellipse-based  DEM  for  investigating  the  uifluence  of  fabric 
anistropy  on  the  mechanical  behavior  of  granular  systems.  Numerous  particulate  assemblages  were  generated 
with  variations  in  particle  angularity  (aspect  ratio),  preferred  bedding  orientation  and  initial  void  ratios.  The 
samples  were  isotropically  compressed  and  then  biaxially  sheared  to  approximately  10%  axial  strain.  This 
Sectiai  describes  the  system  genoaticxi,  fixmation  and  test  simulation  procedures  used  in  the  current  DEM  study. 
The  numerical  samples  were  generated  using  the  preprocessor  program  FEGEN.  Geometry  parameters  and 
matoial  pr(^)erties  were  selected  to  obtain  particulate  assemblages  similar  to  those  used  by  Konishi  et  al  ( 1983). 
Details  of  the  sample  preparation  and  test  results  are  also  foimd  in  Meachum  (1994). 

6.1  Sample  Generation 

Numoous  sanq)les  woe  ^nerated  consisting  of  qjproximately  700  particles  having  aspect  ratios  a/b  (major 
semi-axis  to  mincx  soni-axis)  of  1.0, 1.1, 1.4, 1.7, 2.0, 2.5  and  3.0:1  and  preferred  orientations  0  (measured  from 
the  hOTizontal)  of  0,  30,  60  and  90°.  To  maintain  a  preferred  bedding  plane  during  the  sample's  formation, 
particle  bedding  angles  0  were  assigned  a  standard  deviation  of  ±15°  within  which  the  orientations  could  vary 
sli^tly  while  maintaining  overall  particle  alignment.  At  least  two  samples  were  generated  for  each  combination 
of  a/b  and  0.  Larger  samples  consisting  of  approximately  2000  and  5000  particles  were  also  generated.  Some 
samples  were  generated  without  preferred  particle  alignment  producing  geometric  fabrics  with  more  random 
distributions  of  particle  orientations.  The  gradation  for  all  of  the  numerical  assemblages  consisted  of  an  equal 
number  of  three  different-sized  particles,  the  largest  approximately  two  times  as  large  as  the  smallest  particle. 
Table  6  presents  the  particle  dimensions  of  the  three  different  particle  sizes  for  each  aspect  ratio  group. 

The  geometry  of  the  two-dimensicMial  numerical  assemblages  consisting  of  approximately  700  particles  was 
similar  to  that  us^  Konishi  et  al. ,  measuring  approximately  0.2 1  m  wide  and  0.33  m  high.  The  larger  2000- 
and  5000-particle  samples  were  scaled  proportionally  so  that  the  sample  height  was  approximately  1 .5  times  the 
width. 

Several  packing  schemes  were  used  to  form  numerical  samples  with  variations  in  the  initial  void  ratio. 
Forming  samples  with  loose  packing  proved  very  difficult  as  the  curvature  of  the  ellipses  tended  to  cause  the 
sample  to  be  unstable,  collapsing  into  a  dense  packing  even  under  the  application  of  low  confining  stresses. 

•  O  Sample  Series  -  Particles  were  individually  placed  by  hand  to  obtain  densely  packed  samples  yielding 
void  ratios  similar  to  those  reported  by  Konishi.  The  larger  samples  consisting  of 2004  and  5115  particles 
were  formed  as  composites  of  the  smaller  700-particle  samples. 

•  RSample  Series  -  Scxne  samples  woe  fexmed  by  generating  the  particles  in  a  much  larger  "box"  and  simply 
moving  die  top  and  right-side  rigid  boundaries  (or  walls)  in  using  confining  stresses  of  0.7  kPa.  Control  of 
the  preferred  particle  bedding  to  within  a  standard  deviation  of  ±15°  was  lost,  however,  due  to  particle 
rolling.  The  renting  cxientation  fabrics  were  more  isotropic  than  the  extremely  anisotropic  bedding  planes 
of  the  hand-packed  O  sample  series.  Looser  packing  was  obtained  using  this  packing  scheme  as  more  edge- 
to-face  particle  contacts  existed,  particularly  in  samples  containing  the  flatter  particles. 

•  L  Sample  Series  -  To  form  loosely  packed  numerical  samples  while  maintaining  preferred  horizontal 
particle  bedding  (0  =  0°),  the  particles  were  individually  and  loosely  placed  into  a  sample  box  of  fixed  width 
of  0.21  m  and  varied  height.  A  uniaxial  stress  of  0.7  kPa  was  then  applied  vertically  to  the  assemblage 
producing  a  rather  stable  sample  confined  laterally  by  the  left  and  right  rigid  walls.  After  simulating  for  1 
second,  the  assemblages  were  isotropically  compressed  at  a  mean  stress p  of  0.7  kPa. 
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•  FL  FD  Sample  Series  -  Another  technique  used  to  form  loose  samples  while  maintaming  preferred 
particle  bedding  was  to  vaiy  the  interparticle  sliding  friction  angle  during  consolidation.  Previously  hand- 
packed  Oval  3  and  Oval  5  (a/b  =  2.0)  samples  were  given  a  value  for  <t>^  of  5°  or  52®  during  the  isotropic 
compression  simulation.  When  the  samples  reached  static  equilibrium,  <|>j,  was  returned  to  26°  and  the 
samples  were  then  isotropically  compressed  frirther  at  p  =  0.7.  The  FL  sample  series  are  looser  samples 
using  <|>j,  =  52°  during  formation,  while  the  FD  sample  series  are  denser  samples  using  ^  =  5°  during 
formation. 


Cross-section 

Group 

Aspect  Ratio 
(a/b) 

Particle  Size 

Particle  Diameters  (mm) 
Max/Min 

large 

14.0/14.0 

OvalO 

1.0 

medium 

10.0/10.0 

small 

6.0/6.0 

large 

14.8/13.4 

Ovall 

1.1 

medium 

9.9/8.9 

small 

6.3/5.7 

large 

16.0/11.3 

Oval  2 

1.4 

medium 

10.7/7.4 

small 

7. 1/4.9 

large 

19.52/9.76 

Oval  3 

2.0 

medium 

13.04/6.52 

small 

8.7/4.34 

large 

17.0/10.0 

Oval  4 

1.7 

medium 

12.0/7.0 

small 

8.0/4.8 

large 

18.4/9.2 

Oval  5 

2.0 

medium 

12.8/6.4 

small 

8.8/4.4 

large 

20.4/8.2 

Oval  6 

2.5 

medium 

14.4/5.8 

small 

10.0/4.0 

large 

22.8/7.6 

Oval? 

3.0 

medium 

15.6/5.2 

small 

10.8/3.6 

Table  6.  Dimensions  of  DEM  particles 


52 


Material  Properties  for  DEM  Particles 

Mass  Densi^  p  (g/cm^) 

Poissoii  Ratio  v 
CoefHdent  of  Restitutim  e 
Cdiesion  C  OsN/m^ 

N^odnal  Spring  Stilus  (kN/m^) 

Tan^itial  Spring  Stiffiiess  k.  (kN/m^ 

Sliding  Friction  Angie  O 

Material  Properties  for  DEM  Boundaries  (Walls) 


Coefficient  of  Restitution  e  0,50 

Ck^sicm  C  (kH/m^)  0 

Ncamal  Spring  Stiffness  k„  (kN/m^  30000 

Tai^eotiid  Spring  Stiffiiess  k,(kNr/m^  20000 

Sliding  Fricticai  Angle  (|)^  (*)  0 


1.50 

0.25 

0.50 

0 

2000 
26  or  52 


Table  7.  Material  properties  used  for  the  DEM  simulations. 


All  of  the  numerical  samples  were  generated  having  the  same  material  and  system  prq>erties  and 
gradaticm  while  the  particle  shapes  and  orioitaticMis  and  the  sample  packing  varied.  Table?  presents  the  material 
property  values  used  to  define  the  numerical  systems  in  the  preprocessor  program  FEGEN.  To  eliminate 
unbalanced  fixces  due  to  the  sample's  self-weight,  the  simulations  were  performed  without  gravity. 


6.2  Isotropic  Compression  Simulation 

After  the  samples  wae  formed,  they  were  isotropically  compressed  to  a  mean  stress  />  of  0.7  kPa  using 
the  simulation  program  PARSE.  Higher  isotropic  confining  stresses  of  1 .4  kPa  and  2. 1  kPa  were  also  applied 
to  Oval  2  (a^  =  1.4)  and  Oval  5  (a/b  =  2.0)  samples  with  preferred  horizontal  bedding  (0  =  0°).  Figure  40 
illustrates  a  typical  numerical  sample  isotropically  compressed  at  a  mean  confining  stress  p  of  0.7  kPa.  The 
resulting  normal  contact  force  chains,  shown  in  Figure  4 1,  are  plotted  fi'om  particle  center  to  particle  center  with 
the  relative  thickness  of  the  force  chains  proportional  to  the  magnitudes  of  the  normal  forces. 

The  isotrc^ic  compression  tests  were  simulated  for  1  secoid  until  the  samples  achieved  static  equilibrium 
at  the  desired  mean  confining  stress p  of  0.7  kPa.  The  simulation  time  necessary  to  reach  equilibrium  conditions 
can  vary  depending  cm  the  num^r  of  particles  in  the  samples.  The  smaller  700-particle  sample  achieved 
equilibrium  faster  than  the  larger  5 1 15-particle  sample.  For  the  700-particle  samples,  a  1-second  simulation 
translated  to  approximately  30  hours  of  dedicated  computer  time  on  the  a  DEC  5000/120  UNIX  workstation. 
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02111 
a/b  =  1.4 
0  =  0" 


Figure  40.  Typical  isotropically  compressed  DEM  sample  with  a/b  —1.4  and  0  —  0°  prior  to  shear. 


-02111 
a^  =  1.4 
0  =  0° 


Figure  41.  Face  diains  of  a  typical  isotropically  compressed  DEM  sample  with  a/b  —  1.4  and  0  —  0°  prior  to 
shear 
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6.3  Biaxial  Compression  Test  Simulation 

The  isotropically  compressed  samples  were  sheared  b^ond  failure  under  biaxial  compression  to 
approximately  10%  axial  strain  using  the  simulation  program  PARSE.  The  compression  tests  were  simulated 
moving  the  top  and  bott(xn  walls  toward  each  other  at  a  constant  velocity  ^^le  die  right  side  wall  maintained 
a  cmifining  stress  of  0.7kPa.  Figures  42  and  43  show  the  test  sample  from  Figures  40  and  41  at  approximately 
10%  axial  strain. 

Owing  to  a  need  to  conserve  computational  effort,  it  was  necessary  to  minimize  the  number  of  timesteps 
needed  to  oxiqilete  the  shearing  of  the  sample.  This  entailed  shearing  the  samples  at  the  fastest  possible  strain 
rate,  with  the  largest  possible  time  step  of  integration.  If  the  walls  are  moved  too  fast  (too  hi^  a  strain  rate), 
excessive  ^namic  and  possibly  abnormal  niunerical  effects  may  occur,  such  as  excessive  overlapping  between 
neighboring  particles.  It  is  desirable  to  move  the  walls  as  fast  as  possible  without  significantly  compromising 
the  model's  accuracy. 

To  determine  how  fast  the  biaxial  shear  could  be  conducted  without  compromising  accuracy,  four  biaxial 
ccxnpressicn  tests  were  pofcxmed  on  the  same  sanqile  with  diffaent  wall  velocities.  To  achieve  10%  axial  strain 
fcx  each  test  simulation,  top  (-)  and  bottom  (+)  wdl  velocities  of  ±0.0192,  ±0.0096,  ±0.0064  and  ±0.0048  m/s 
were  used  to  compress  the  sanqiles  for  1, 2, 3  and  4  seconds,  respectively.  Figures  44  and  45  plot  the  deviator 
stress>axial  strain  and  volumetric  strain-axial  strain  data  for  these  tests,  which  indicate  little  difference  in  the 
observed  overall  strength  and  deformation  behavior.  Based  on  these  results,  the  numerical  assemblages  were 
biaxially  compressed  using  the  faster  simulation  time  of  1  second  with  wall  velocities  of  ±0.016  m/s,  which 
iqjresents  an  axial  strain  rate  of  about  0. 1  s'* .  Note  that  in  Figures  44,  the  deviator  stress  is  quantified  using  the 
n^ilized  friction  angle  which  was  previously  defined  in  equation  (56).  Also,  note  that  in  Figure  45  and  in  all 
subsequent  plots,  volumetric  strain  is  plotted  positive  for  contraction  and  negative  for  dilation  (even  though  this 
is  opposite  to  the  compression  positive  sign  convention). 

To  detomine  how  large  a  timestep  can  be  used,  several  considerations  need  to  be  made.  The  ellipse-based 
DEM  model  suggests  a  timestep  At  based  on  several  considerations  including  the  natural  vibration  of  a  single 
spring-mass  ^rston,  multiple  spring-single  mass  system,  and  ensuring  that  colliding  particles  remain  in  contact 
fcx  at  least  thiw  timestqis.  Th^  tests  are  described  in  d^ain  in  (^<xkum  and  Ting  ( 1 986).  For  the  given  particle 
geaneby  and  contact  properties  used  in  the  current  study,  a  timestep  of  2  x  10'^  s  was  suggested  by  FEGEN  and 
PARSE.  To  test  wheto  the  suggested  timestep  was  adequate  to  produce  reasonably  accurate  results,  a  numerical 
sample  was  biaxially  sheared  under  the  same  simulation  conditions  except  with  variations  in  the  timestep.  The 
results  for  the  nusbilized  internal  friction  angle  ^  and  the  volumetric  strain  of  these  simulations  are  shown  in 
Figures  46  and  47.  The  minor  behavioral  differences  using  the  different  timesteps  suggests  that  use  of  At  =  2 
X  10'^  s  is  appropriate. 
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Figure  42.  Biaxially  compressed  DEM  sample  with  a/b  =  1.4  and  0  =  0®  at  approximately  10%  axial  strain 


Figure  43.  Fcxice  chains  fix  biaxially  ccxnpressed  DEM  san^le  with  a/b  =1.4  and  0  —  0°  at  approximately  1 0% 
axial  strain 
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1 -second  simulation 

2-second  simuiation 

3-second  simulation 

-w-  4-second  simuiation 

Figure  44.  Mobilized  internal  friction  angle  for  a/b  =  3.0  and  0  =  0°  sample  biaxially  sheared  using  different 
wall  velocities 


1 -second  simulation 

2-second  simulation 

3-second  simulation 

4-second  simulation 

Figure  45.  Volumetric  strain  behavior  for  a/b  =  3.0  and  0  =  0°  sample  biaxially  sheared  using  different  wall 
velocities 


57 


0  2  4  6  8  10 


Axial  Strain  (%) 


Figure  46.  Mobilized  internal  angle  of  friction  fw  a/b  =  3.0  and  0  =  0°  sample  biaxially  sheared  using  different 
timesteps  dt. 
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Figure  47.  Volumetric  strain  behavior  for  a/b  =  3.0  and  0  =  0°  sample  biaxially  sheared  using  different 
timesteps  dt. 
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7.  COMPARISON  WITH  KONISHI,  ODA  AND  NEMAT-NASSER  RESULTS 


As  previously  described  in  Section  5.4,  Konishi,  Oda  and  Nemat-Nasser  (1983,  1982)  investigated  the 
influence  of  particle  shape  on  the  overall  mechanical  behavior  of  granular  materials  using  photoelastic 
polyurethane  rubber  rods.  In  that  Section,  a  series  of  detailed  photographic  plates  supplied  by  Prof  Oda  was 
digitized  and  used  as  itq)ut  to  the  DEM  simulaticxi.  This  single  ph^ical  test  was  then  used  in  the  validation  phase 
of  this  stu^.  In  the  curroit  Secdai,  the  other  expoimaital  results  from  the  Konishi  et  al  tests  are  compared  with 
DEM  results.  Hie  numerical  tests  in  this  Section  have  the  same  overall  general  geometry  and  particle  gradation 
as  the  physical  tests,  but  not  the  exact  initial  locations.  As  before,  the  DEM  particles  have  slightly  different 
shapes  (ellipse  vs  oval)  but  the  same  overall  dimensions.  Results  obtained  from  the  ellipse-based  DEM  test 
simulations  are  similar  qualitatively  to  the  physical  test  data  reported  by  Konishi  et  al.  The  numerical  results, 
howevCT,  yield  consistently  lower  stresses  than  the  physical  test  results.  As  previously  described  for  the  single 
test  conqiarison  in  Secticxi  5.6,  it  is  coijectured  that  these  differences  are  due  to  differences  in  the  particle  shapes. 


7.1  Data  Comparison  for  Samples  with  <|>^  =  26° 

Figure  48  shows  the  influence  of  particle  angularity  and  preferred  (or  aligned)  particle  bedding  on  the 
peak  internal  frktkxi  angle  <!>„„„(  for  the  physical  and  DEM  samples  with  interparticle  sliding  friction  angle  i|>^  of 
26°.  Qualitatively,  the  results  are  similar.  In  general,  higher  values  for  <1>„^  are  achieved  by  samples  with 
preferred  hcaizcmtal  bedding  (0  =  0°),  while  fix  samples  with  60  and  90°  bedding,  the  values  are  lower.  Samples 
consisting  of  angular  particles  (a/b  =  1.4)  yield  higher  internal  friction  angles  than  those  consisting  of  more 
rounded  particles  (a/b  =  1.1).  The  results  of  the  physical  tests  performed  by  Konishi  et  al,  however,  are 
significantly  higher  than  the  results  obtained  from  the  DEM  test  simulations. 


-  DEM  data  (a/b  s  1.1)  a  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  48.  Comparison  of  DEM  results  and  Konishi  et  al  (1983)  physical  results  for  samples  with  an 
interparticle  sliding  friction  angle  (t>^  of  26° 
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The  peak  frictiai  angles  and  initial  void  ratios  (e^)  for  the  physical  and  numerical  samples  are  presented 
in  Table  8.  The  DEM  results  are  average  values  fw  numerical  samples  having  the  same  particle  aspect  ratio  a/b, 
bedding  plane  0  and  confining  stress  (oj  =  0.7  kPa).  In  general,  jL  increases  with  particle  an^larity,  the 
excepticxi  being  the  physical  sanqile  with  60®  bedding.  The  slightly  looser  numerical  samples  yield  significantly 

lower  values  for  ()>„,„  than  the  physical  samples  used  by  Konishi. 
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Table  8.  Results  of  biaxial  compression  tests  performed  by  Konishi  et  al  and  DEM  test  simulations 


Figures  49  to  52  show  the  progressiwi  of  the  principal  stress  ratio  CT,/a2  during  the  course  of  deformation  for 
the  pltysical  and  numerical  samples  with  bedding  planes  of  0, 30, 60  and  90®,  respectively.  The  data  yielded  by 
both  the  physical  and  numerical  samples  show  that  a,/c2  achieves  higher  peak  values  for  samples  with  angular 
particles  (a/b  =  1.4).  The  peak  principal  stress  ratio  achieved  by  the  physical  samples  is  significantly  higher  than 
the  DEM  samples  and  is  achieved  at  much  lower  failure  strains.  The  differences  in  the  residual  stress  ratios, 
howevCT,  are  rather  insignificant.  The  dffieraices  in  peak  stress  ratio  between  the  Konishi  data  and  the  DEM  data 
are  worst  in  samples  with  horizontal  bedding  (0  =  0®). 
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-  OEM  data  (a/b=  1.1)  a  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =1.1)  -a-  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  49.  Con^ariscm  of  DEM  and  Konishi  et  al  physical  results  for  samples  with  0°  bedding  and  inteiparticle 
sliding  friction  angle  (|>^  of  26°. 


Figure  50.  Comparison  of  DEM  and  Konishi  et  al  physical  results  for  samples  with  30°  bedding  and 
interparticle  sliding  friction  angle  <|>^  of  26°. 
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-  DEM  data  (a/b=  1.1)  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =  1 .1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  51.  Comparison  of  DEM  and  Konishi  et  al  physical  results  for  samples  with  60°  bedding  and 
interparticle  sliding  friction  angle  4)^  of  26°. 
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DEM  data  (a/b  =  1.1) 

Konishi  et  al  data  (a/b  =  1.1) 


DEM  data  (a/b  =  1.4) 

Konishi  et  al  data  (a/b  =  1 .4) 


Figure  52.  Comparison  of  DEM  and  Konishi  et  al  physical  results  for  samples  with  90°  bedding  and 
interparticle  sliding  friction  angle  <|>^  of  26°. 
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Figures  53  to  56  plot  the  physical  and  DEM  volumetric  strain  AVA^^  data  for  assemblages  with  0, 30, 60 
and  90°  bedding  (with  dilation  plotted  positive).  Assemblages  with  angular  particles  (a/b  =  1.4)  and  preferred 
horizontal  bedding  (0  =  0°)  initially  ctmtract  mwe  and  longer  than  those  consisting  of  m(M«  roimded  particles 
(a/b  =  1 . 1),  but  ultimately  dilate  more  and  faster  as  deformation  continues.  For  assemblages  with  60°  bedding, 
however,  particle  angularity  has  the  qjposite  elOfect  cm  AV/V^.  Samples  consisting  of  angular  particles  experience 
less  dilation  throughout  the  course  of  deformation  than  those  wi&  rounded  particles.  Though  the  physical  and 
DEM  sanqrles  both  exhibit  volumetric  strain  behavior  similar  qualitatively  to  that  of  real  granular  materials,  the 
physical  samples  are  stiffer,  dilating  significant^  more  than  the  numerical  samples  and  reaching  their  peak  values 
at  much  lower  axial  strains.  The  differences  in  the  Konishi  and  DEM  results  for  AV/V^  are  larger  for  samples 
with  horizontal  bedding,  while  for  samples  with  60°  bedding,  the  differences  are  smaller. 


■  DEM  data  (a/b  s  1.1)  a  DEM  data  (a/b  -  1 .4) 

-B-  Konishi  et  al  data  (a/b  =1.1)  Konishi  et  al  data  (a/b  =  1 .4) 


Figure  53.  Ccmqjariscm  of  DEM  and  Kcmishi  et  al  physical  results  fix’  samples  with  0°  bedding  and  interparticle 
sliding  fiiction  angle  <|>^  of  26°. 
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■  DEM  data  (a/b=  1.1)  ^  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  54.  Comparison  of  DEM  and  Konishi  et  al  physical  results  for  samples  with  30°  bedding  and 
interparticle  sliding  friction  angle  of  26°. 


-  DEM  data  (a/b  =  1.1)  ^  DEM  data  (a/b  =  1 .4) 

Konishi  et  al  data  (a/b  =1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  55.  Comparison  of  DEM  and  Konishi  et  al  physical  results  for  samples  with  60°  bedding  and 
interparticle  sliding  friction  angle  <|>ji  of  26°. 
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■  DEM  data  (a/b=  1.1)  a  DEM  data  (a/b=  1.4) 

-a-  Konishi  et  al  data  (a/b  =  1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  56.  Comparison  of  DEM  and  Konishi  et  al  physical  results  for  samples  with  90°  bedding  and 
interparticle  sliding  frictimi  angle  <|>^  of  26° 


7.2  Data  Comparison  for  Samples  with  =  52° 

Figure  57  plots  the  influence  of  particle  angularity  and  fneferred  particle  bedding  on  the  peak  internal  friction 
angle  <t>™x  for  the  physical  and  DEM  samples  with  interparticle  sliding  friction  angle  <|)i  of  52°.  Higher  values 
for  <|»^  are  achieved  by  samples  with  preferred  horizontal  bedding  (0  =  0°),  while  for  samples  with  60  and  90° 
bedding,  (|>„,„  is  lower.  Assemblages  with  angular  particles  (a/b  =  1 .4)  yield  higher  peak  internal  friction  angles 
than  those  with  more  rounded  particles  (a/b  =  1.1).  As  with  the  results  for  <|)j,  =  26°,  the  physical  results  are 
qualitatively  similar  to  the  DEM  results,  but  quantitatively  much  larger  in  strength. 

Table  9  presoits  the  peak  friction  angle  (|>n„^  and  initial  void  ratio  Oq  for  the  physical  and  numerical  samples 
with  <|)^ = 52°.  The  DEM  results  are  average  values  for  numerical  samples  having  the  same  particle  aspect  ratio 
a/b,  bedding  plane  0  and  confining  stress  (02  =  0.7  kPa).  Both  the  Konishi  et  al  and  DEM  assemblages  with  <|>^ 
=  52°  exhibit  higb^  stiaigth  than  those  with  (|>^  of  26°.  As  with  the  results  for  =  26°,  <|>„„  increases  as  particle 
angularity  increases,  the  exception  being  the  numerical  sample  with  90°  bedding.  The  differences  between  the 
physical  and  numerical  results,  however,  are  still  very  significant.  Though  the  relative  densities  of  the  physical 
and  DEM  san^les  are  similar,  the  test  simulations  yield  significantly  lower  values  for  than  the  physical  tests. 
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■  DEM  data  (a/b=  1.1)  ^  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =  1 .1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  57,  C<»q>arison  of  DEM  and  Konishi  et  al  pltysical  results  fw  samples  with  interparticle  sliding  friction 
angle  <t»^  of  52°. 
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♦  Average  values  for  numerical  samples  with  similar  geometries 

Table  9.  Results  of  biaxial  compression  tests  performed  by  Konishi  et  al  and  DEM  test  simulations. 
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Figures  58  to  61  plot  the  principal  stress  ratio  o,/o2  vs  axial  strain  for  the  physical  and  DEM  tests  with  ()»^ 
of  52®  and  bedding  planes  of  0, 30, 60  and  90®,  respectively.  Though  the  peak  values  achieved  by  the  physical 
and  numerical  tests  do  not  agree,  the  residual  values  compare  favorably.  Also,  the  residual  values  for  o,/o2  are 
about  the  same  regardless  of  differences  in  particle  angularity  and  orientation.  As  with  <|>^  =  26®,  the  physical 
samples  with  <|>j, = 52°  exhibit  stiffer  behavior,  yielding  higher  peak  values  for  C1/02  at  lower  failure  strains  than 
the  DEM  samples. 


•  DEM  data  (a/b  =  1 .1)  ^  DEM  data  (a/b  =  1 .4) 

Konishi  et  al  data  (a/b  =1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  58.  Conq)arison  of  DEM  and  Ktmishi  et  al  results  for  samples  with  0®  bedding  and  interparticle  sliding 
fiiction  angle  <j)^  of  52°. 
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-  DEM  data  (a/b  =  1.1)  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =  1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  59.  Comparisrxi  of  DEM  and  Kcxiishi  et  al  results  for  sanqjles  with  30®  bedding  and  interparticle  sliding 
friction  angle  <|>^  of  52®. 
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-  DEM  data  (a/b  =  1.1)  a  DEM  data  (a/b  =  1.4) 

-B-  Konishi  et  al  data  (a/b  =  1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  60.  Qxnparison  of  DEM  and  Konishi  et  al  results  for  sanples  with  60®  bedding  and  interpaiticle  sliding 
friction  angle  <1>^  of  52®. 
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■  DEM  data  (a/b=  1.1)  a  DEM  data  (a/b=  1.4) 

Konishi  et  al  data  (a/b  =  1.1)  Konishi  et  al  data  (a/b  =1.4) 

Figure  61.  Cotnqparistxi  of  DEM  and  Kcmishi  et  al  results  frx  samples  with  90*^  bedding  and  interparticle  sliding 
friction  angle  of  52°. 


Figures  62  to  65  plot  the  volumetric  strain  behavior  for  the  physical  and  numerical  tests  with  an  interparticle 
sliding  friction  angle  of  52°  and  0,  30,  60  and  90°  bedding,  respectively.  Samples  consisting  of  angular 
particles  (a/b  =  1.4)  with  preferred  horizontal  bedding  (0  =  0°)  experience  more  and  Icmger  initial  contraction 
than  diose  with  more  roun^  particles  (a^  =1.1).  At  larger  axial  strains,  however,  assemblages  with  angular 
particles  dilate  mcxe.  Both  die  physical  and  DEM  assemblages  with  a/b  =  1.4  and  60°  bedding  dilate  less  during 
the  course  of  deformation  than  Aose  with  a/b  =  1. 1.  Though  both  the  physical  and  numerical  samples  exhibit 
volumetric  strain  behavior  similar  qualitatively  to  that  of  real  granular  materials,  the  physical  samples  are  stiffer, 
dilating  signifrcantly  more  and  faster  than  the  DEM  samples. 
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■  DEM  data  (a/b  =  1.1)  ^  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =  1 .1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  62.  Comparison  of  DEM  and  Kcxiishi  et  al  results  for  samples  with  0°  bedding  and  interparticle  sliding 
friction  angle  <t>n  of  52®. 


-  DEM  data  (a/b  =  1 .1)  ^  DEM  data  (a/b  =  1 .4) 

Konishi  et  al  data  (a/b  =1.1)  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  63.  Comparistxi  of  DEM  and  Kcmishi  et  al  results  fix  samples  with  30®  bedding  and  interparticle  sliding 
firiction  angle  <))j.  of  52®. 
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-  DEM  data  (a/b  =  1.1) 

A  DEM  data  (a/b  =  1.4) 

-B-  Konishi  et  al  data  (a/b  =  1.1) 

-A-  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  64.  Compariscm  of  DEM  and  Konishi  et  al  results  for  samples  with  60°  bedding  and  interparticle  sliding 
friction  angle  <t>^  of  52°. 


-  DEM  data  (a/b  =1.1)  ^  DEM  data  (a/b  =  1 .4) 

-B-  Konishi  et  al  data  (a/b  =1.1)  -a-  Konishi  et  al  data  (a/b  =  1 .4) 

Figure  65.  Ccxnparison  of  DEM  and  Kraiishi  et  al  results  for  samples  with  90°  bedding  and  interparticle  sliding 
friction  angle  of  52°. 
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7J  Discussion  of  Data  Comparison 

As  with  the  validaticxi  test  on  the  single  digitized  sample,  comparison  of  the  physical  and  numerical  biaxial 
shear  test  results  using  san^les  with  similar  gradation  and  overall  properties  demonstrated  significant  differences 
in  strength  and  deformation  behavica'.  Hi^er  strength  and  more  dilation  are  achieved  by  the  physical  samples 
at  much  lower  axial  strains.  The  discrepancies  in  the  physical  and  DEM  data  may  be  caused  by  a  number  of 
reasons.  The  accuracy  of  the  ellipse>based  DEM  model,  the  nature  of  the  contact  interaction,  the  selection  of 
properties  and  parameters  and  the  differences  in  the  shape  of  the  particles  (oval  V5.  ellipse)  are  considered  as 
possible  causes  for  the  discrepancy  in  results. 

Based  on  extensive  testing  of  the  code  as  outlined  in  Section  4.4,  one  presumes  that  the  ellipse-based  DEM 
model  is  working  ]H<^)erly.  Based  on  die  results  described  in  Section  5.5,  variation  of  the  form  and  value  of  the 
contact  properties  doM  not  ^pear  to  significantly  alter  the  observed  macroscopic  strength  and  volume  behavior. 
An  alternate  explanation  for  the  discrepancy  is  based  on  subtle  differences  in  the  shape  of  the  physical  and 
numerical  particles.  This  hypothesis  was  examined  for  the  single  digitized  test  in  Section  5.6.  In  this  Section, 
this  hypothesis  is  further  examined  for  the  tests  on  samples  at  different  bedding  angles. 

As  mentioned  earlier,  Kcxiishi  et  al  used  photoelastic  rubber  rods  with  oval-shaped  cross-sections,  while  the 
current  study  uses  ellipse-shaped  particles.  A  close  visual  examination  of  the  oval-shaped  rods  in  Figures  29 
shows  that  their  sides  are  flattened,  while  the  surface  of  the  ellipse-shaped  particles  are  continuously  rounded. 
Presumably,  the  flattened  sides  facilitate  a  more  stable  packing  when  the  particles  are  assembled  or  "stacked". 
Whfn  subjected  to  a  vertical  shear  stress,  the  physical  samples  form  particulate  columns  that  are  far  more  stable 
and  can  resist  significantly  higher  shear  stresses  than  the  ellipse-based  assemblages,  especially  when  the  bedding 
plane  is  hori7ontal  Foe  the  1.4: 1  aspect  ratio  test,  this  is  analogous  to  axial  loading  of  an  imcemented  brick  wail. 
The  particulate  columns  formed  by  the  ellipse-shaped  particles,  however,  are  less  stable  and  hence  weaker. 

Figures  66  and  67  compare  the  rose  diagrams  of  the  contact  normal  distribution  prior  to  shear  and  at  peak 
strength  for  the  numerical  and  physical  samples  with  a  particle  aspect  ratio  of  1.4,  preferred  horizontal 
bedding  (0  =  0°)  and  interparticle  sliding  fhetion  angle  (|>^  of  26°.  The  contact  statistics  arc  computed  assuming 
two  cmtact  normals  exist  at  eadi  ccxitact  The  rosettes  for  the  numerical  assemblages  consider  only  interparticle 
contacts  and  exclude  wall-particle  contacts.  The  physical  sample  yields  a  more  amsotropic  distribution  prior  to 
shear  achieving  a  maximum  fiequetx^  of  ^proximately  13%,  while  the  numerical  sample  is  less  anisotropic  and 
adiieves  a  maximum  fiequen^^  of  oily  9%.  This  suggests  that  the  physical  sample  is  mOTe  stable  prior  to  shear 
as  more  contact  ncamals  are  oriented  toward  the  vertical  and  arc  predisposed  to  provide  more  resistance  to  the 
vertically  tqq)lied  shear  stress  Oj.  At  peak  strength,  the  contact  normal  distribution  for  Konishi's  assemblage  is 
again  more  anisotrc^ic  than  the  numerical  sample  as  the  flattoied  sides  of  the  oval-shaped  particles  produce  more 
contact  ncHrmals  that  are  coaxial  with  a,. 

Disorqjancies  also  exist  in  the  data  comparison  for  samples  having  preferred  bedding  planes  of  30, 60  and 
90°.  These  differences,  however,  are  not  as  large  as  those  observed  for  samples  with  horizontal  bedding.  The 
smallest  difference  in  occur  in  samples  with  60°  bedding  with  Konishi  et  al  reporting  a  value  of  39°  and  the 
DEM  yielding  an  average  value  of  30°.  Presumably,  the  flattened  sides  of  the  rubber  particles  become  less 
influential  on  tiie  Sample's  stroigth  behavicx'  as  the  orientaticxi  of  the  contact  plane  and  the  direction  of  the  applied 
shear  stress  beccanes  skewed. 

Figures  68  and  69  compares  the  contact  normal  fl'cquency  distributions  prior  to  shear  and  at  peak  strength 
for  physical  and  DEM  samples  with  a^  =  1.4, 0  =  60°  and  <|>^  =  26°.  While  the  rosettes  show  that  the  oval- 
shaped  particles  possess  a  more  skewed  distribution  of  contact  normals,,  the  differences  are  less  sigmficant  than 
those  observed  in  samples  with  preferred  horizontal  bedding  (0  =  0°). 
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Figure  66.  FrequeiK^  distributicm  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -  1.4, 
0  =  0®  and  <t)^  =  26®  prior  to  shear. 


Figure  67.  Frequency  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1.4, 
0  =  0®  and  <t>^  =  26®  at  peak  strength. 
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Similar  results  are  obsoved  in  axipariscais  of  the  contact  normal  distributions  for  the  remaining  Oval  1  and 
Oval  2  sanqjles.  Figures  70  to  77  show  ccunparisons  of  contact  normal  distributions  for  the  Oval  1  (a/b  =1.1) 
samples  with  preferred  bedding  planes  of 0,  30, 60  and  90°,  respectively,  and  an  interparticle  sliding  friction 
angle  of  26°.  Figures  78  to  81  show  the  same  comparison  for  the  Oval  2  (a^  =  1.4)  samples  with  30  and  60° 
bedding,  respectively.  Gxtq>arisons  of  the  contact  normal  distributions  for  samples  with  (j>^  of  S2°  are  shown  in 
Figures  82  to  97.  Based  on  the  distribution  of  the  contact  normals,  it  is  likely  that  the  discrepancies  observed  in 
the  strength  behavior  of  the  physical  and  numerical  samples  are  due  to  the  dissimilarities  in  the  particle  shapes. 


Figure  70.  Frequency  distributi<xi  of  ctmtact  ncmnals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 . 1, 
0  =  0°  and  <|>^  =  26°  prior  to  shear. 


Figure  71 .  Frequoxy  distributicm  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 . 1 , 
0  =  0°  and  <|>^  =  26°  near  peak  strength. 
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Figure  72.  Frequracy  distributicxi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -1.1 
0  =  30®  and  (j)^  =  26®  prior  to  shear. 
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Figure  73.  Frequaicy  distributicMi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -  1 . 1 
6  =  30®  and  <1>^  =  26°  near  peak  strength. 


Figure  74.  Frequency  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 . 1 , 
0  =  60®  and  (t>^  =  26®  prior  to  shear. 


Figure  75.  Frequency  distributicm  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 . 1 , 
0  =  60®  and  <t)^  =  26®  near  peak  strength. 


Figure  76.  Frequaxy  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a^  =  1. 1 
0  =  90°  and  (}>,,  =  26°  prior  to  shear. 


Figure  77.  Frequoicy  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  1 . 1 
0  =  90°  and  (t>„  =  26°  near  peak  strength. 


Figure  78.  Frequency  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 .4, 
0  =  30®  and  <t>^  =  26®  prior  to  shear. 


Figure  79.  Frequency  distributicai  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 .4, 
0  =  30®  and  (|)^  =  26®  near  peak  strength. 
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Figure  80.  Frequent  distributiai  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -  1 .4, 
0  =  90°  and  <|>^  =  26°  prior  to  shear. 


Figure  81.  Frequency  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -  1 .4, 
0  =  90°  and  <t»^  =  26°  near  peak  strength. 


Figure  82.  Frequency  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 . 1 , 
0  =  0®  and  <|)^  =  52®  prior  to  shear. 


Figure  83.  Frequency  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 . 1 , 
0  =  0®  and  <1)^  =  52®  near  peak  strength. 
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Figure  84.  Frequency  distributi(xi  of  contact  ncmnals  for  numerical  (left)  and  physical  samples  with  a/b  -  1 . 1 , 
0  =  30°  and  <|)^  =  52°  prior  to  shear. 


Figure  85.  Frequoicy  distribution  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -1.1, 
0  =  30°  and  <|»^  =  52°  near  peak  strength. 
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Figure  86.  Frequency  distributicn  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 . 1 , 
0  =  60®  and  (})^  =  52®  prior  to  shear. 


Figure  87.  Frequency  distributiai  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a^  =  1 . 1, 
0  =  60®  and  (|»^  =  52®  at  peak  strength. 
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Figure  89.  Frequoicy  distributicxi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  —1.1, 
0  =  90®  and  <t>n  =  52°  near  peak  strength. 


Figure  90.  Frequency  distributicm  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a^  =  1 .4, 
6  =  0°  and  <|i^  =  52°  prior  to  shear. 


Figure  91.  Frequoxy  distributicxi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 .4, 
0  =  0°  and  <l»jj  =  52°  near  peak  strength. 
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Figure  92.  Frequraicy  distributiOTi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a^  =  1.4, 
0  =  30®  and  i|>^  =  52°  prior  to  shear. 


Figure  93.  Frequency  distributicxi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a^  -  1 .4, 
0  =  30°  and  =  52°  at  peak  strength. 


Figure  94.  Frequency  distributicxi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  =  1 .4, 
0  =  60°  and  <|)^  =  52°  prior  to  shear. 


Figure  95.  Frequoxy  distributicHi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a^  =  1 .4, 
0  =  60°  and  <t>j.  =  52°  near  peak  strength. 


Figure  96.  Frequency  distributicKi  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -  1 .4, 
0  =  90®  and  (t>^  =  52°  prior  to  shear. 


Figure  97.  FrequoKy  distributiai  of  contact  normals  for  numerical  (left)  and  physical  samples  with  a/b  -  1 .4, 
0  =  90®  and  <1>^  =  52®  at  peak  strength. 
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8.  RESULTS  OF  DEM  ANALYSES 


This  Section  presents  the  results  of  the  biaxial  compression  test  simulations  cm  numerous  cohesionless 
granular  assranblages  with  various  combinations  of  particle  angularity  and  preferred  bedding  plane  mientation 
as  well  as  variations  in  the  initial  void  ratio,  packing  and  confining  stress.  The  interparticle  sliding  fiiction  angle 
for  all  of  the  assemblages  presented  is  26**.  Detailed  simulatiai  results  for  each  numerical  sample  are  provided 
in  Appendix  A  of  Meachum  (1994).  The  strength  and  deformation  results  are  discussed  in  Section  9. 

8.1  Influence  of  Particle  Shape  on  Strength 

Particle  shape  has  been  shown  by  numerous  experimental  and  numerical  studies  to  have  a  strong  influence 
on  the  overall  strength  behavior  of  granular  materials.  Results  show  that  even  slight  increases  in  particle 
angularity  can  significantly  increase  strength.  It  has  been  suggested  that  the  increase  in  strength  is  due  to  the  fact 
that  angular  particles  interlock  more  thereby  inhibiting  the  negative  effects  due  to  particle  rolling  (ting  et  al 
1993).  Coisequently,  particles  must  slide  past  each  other  resulting  in  more  fiictional  shear  resistance  and  higher 
ovei^dl  shear  strength.  Angular  particles  possess  orientation  and  may  form  preferred  (or  aligned)  bedding  planes 
producing  dense  assemblages  with  extremely  anisotropic  fabrics.  Experimental  studies  have  shown  that  strength 
is  significantly  reduced  when  the  orientations  of  the  contact  normal  plane  and  the  major  principal  stress  o,  are 
skewed  (Ladd  et  al  1977). 

The  ellipse-based  DEM  results  show  good  quantitative  and  qualitative  agreement  with  the  results  fi’om 
experimental  studies.  Strength  is  shown  to  increase  as  particle  angularity  increases.  The  strength  decreases 
si^iificantly,  however,  when  the  (xientation  of  the  bedding  plane  is  aligned  toward  the  predicted  Mohr-Coulomb 
failure  plane  (4S‘’  +  (|>/2).  Looser  assemblages  are  shown  to  exhibit  less  strength  than  the  denser  assemblages. 
Finally,  assemblages  with  more  isotropic  fabrics  tend  to  be  looser  and  exhibit  less  strength  than  those  with 
anisotropic  fabrics. 

8.1.1  Effect  of  Particle  Aspect  Ratio 

The  effect  of  particle  aspect  ratio  (a/b)  on  the  peak  internal  fiiction  angle  <t)„ax  is  shown  in  Figure  98.  Note 
that  the  lines  are  plotted  through  the  averages  for  two  samples  with  similar  geometries.  A  significant  increase 
in  is  associated  with  increased  particle  “angularity”.  DEM  samples  consisting  of  "flat"  particles  (a/b  =  3.0) 
yield  values  fw  of  iq)proximately  55°  while  the  san:q)les  consisting  of  round  particles  (a/b  =1.0)  yield  values 

of  only  25°.  Even  a  slight  increase  in  particle  aspect  ratio  is  seen  to  significantly  increase  the  strength  of  the 
assemblage  as  samples  with  a/b  =  1. 1  yield  values  for  <{»„,„  of  approximately  30°.  The  results  of  the  DEM  test 
simulaticms  shown  in  Figure  98  arc  similar  qualitatively  and  quantitatively  to  results  of  physical  tests  performed 
(XI  real  granular  matoials  shown  in  Figure  1 .  Physical  samples  consisting  of  more  angular  Olivine  sand  particles 
yield  significantly  higher  values  for  <|>^  than  the  more  roimded  glass  beads. 

8.1.2  Effect  of  Aspect  Ratio  and  Particle  Bedding 

The  munerical  test  results  shown  in  Figure  98  are  for  samples  with  preferred  horizontal  bedding  (0  =  0°). 
Figure  99  shows  that  <!>,„„  increases  as  the  particles  become  more  flat,  but  when  the  bedding  plane  is  inclined 
toward  the  vertical,  (j)^  decreases  significantly.  For  assemblages  consisting  of  flatter  particles  (a/b  >  2.0), 
is  minimum  wlien  the  bedding  plane  is  cxiented  60°  fi-om  the  horizontal  (or  approximately  45°  +  ^/2).  For  these 
samples,  is  approximately  the  same  as  the  interparticle  fnction  angle  (26°).  For  all  of  the  numerical 
samples  wiA  variatkxis  in  particle  aspect  ratio  and  bedding,  however,  the  peak  internal  fnction  angle  is  minimum 
for  samples  consisting  of  round  particles. 
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■  0°  bedding  | 

Figure  98.  Effect  of  particle  aspect  ratio  on  the  peak  friction  angle  <l>nax  for  DEM  samples  with  preferred 
horizontal  particle  bedding  (0  =  0®). 
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Figure  99.  Effect  of  particle  aspect  ratio  <xi  the  peak  fricti<Mi  angle  for  DEM  samples  with  preferred  particle 

bedding  0. 
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The  effect  of  particle  aspect  ratio  and  orientation  on  the  peak  internal  fricticm  angle  is  also  shown  in 
Figure  100.  Again,  note  that  the  lines  are  plotted  through  the  averages  for  two  samples  with  similar  geometries. 
The  maximum  internal  friction  angle  hence  the  strength,  is  highest  for  samples  with  horizontal  bedding  (0 
=  0®)  and  lowest  for  samples  with  60  and  90®  bedding.  The  DEM  results  for  <|)„^  shown  in  Figure  100  are 
qualitatively  similar  to  the  results  of  physical  tests  performed  on  different  sands  shown  in  Figure  101.  <|)„„  for 
tte  numerical  assemblages  with  60  and  90®  bedding,  however,  is  lower  than  for  the  physical  samples.  Ladd  et 
al  (1977)  defines  p  as  the  angle  between  the  contact  normal  plane  and  the  orientation  of  the  applied  principal 
stress  o,  which  is  in  effect  the  same  as  the  angle  6  used  herein  to  define  the  orientation  of  bedding  plane. 

Similar  peak  strength  behavior  for  the  numerical  samples  is  illustrated  by  the  shear  stress  q  results  shown 
in  Figures  102  and  103  and  the  principal  stress  ratio  o,/a2  results  shown  in  Figures  104  and  105.  A  maximum 
shear  stress  of  ^proximately  4  kPa  is  achieved  by  samples  consisting  of  "flat"  particles  (a/b  =  3.0)  with  0® 
bedding,  Miile  those  consisting  of  round  particles  (a^  =  1.0)  yield  values  for  q„ax  only  0.6  kPa.  Similarly,  a 
maximum  value  fw  O1/C2  of  {q)proximately  10  is  achieved  by  samples  with  a/b  =  3.0  and  0  =  0®,  while  those  with 
a^  =  1.0  yield  values  of  only  2.5.  Values  for  q„,„  and  (ci/a2)Lx  are  considerably  lower,  however,  when  the 
bedding  plane  is  oriented  toward  45°  +  <|)/2.  For  samples  with  30, 60  and  90°  bedding,  the  increase  in  strength 
due  to  particle  aspect  ratio  appears  to  be  less  significant  when  a/b  is  greater  than  1 .4. 
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Figure  100.  Effect  of  prefened  particle  bedding  0  on  the  peak  internal  firiction  angle  (t>n^  for  DEM  samples 
with  variations  in  particle  aspect  ratio  a/b. 
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ANGLE,  p  (Degrees) 

Figure  101.  Effect  of  preferred  particle  bedding  p  on  the  internal  friction  angle  for  different  sands  (from 
Laddctfl/ 1977). 


Figure  102.  Effect  of  particle  aspect  ratio  cm  the  peak  shear  stress  for  DEM  samples  with  preferred 
particle  bedding  6. 
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Particle  Orientation  (deg) 
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Figure  103.  Effect  of  preferred  particle  bedding  6  on  the  peak  shear  stress  for  DEM  samples  with 
variations  in  particle  aspect  ratio  a/b. 
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Figure  104.  Effect  of  particle  aspect  ratio  on  the  peak  principal  stress  ratio  (cr,/a2)mn  for  DEM  samples  with 

preferred  particle  bedding  6. 
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Figure  105.  Effect  of  preferred  particle  bedding  0  on  the  peak  principal  stress  ratio  for  DEM 

samples  with  variations  in  particle  aspect  ratio  a^. 


Figures  106  and  107  illustrate  the  effect  of  particle  aspect  ratio  on  the  progression  of  the  mobilized  internal 
friction  angle  <|)  during  the  course  of  deformation  for  samples  with  0  =  0  and  60°,  respectively.  As  shown 
previously,  tire  eUipse-based  DEM  model  ccxnputes  mobilized  <|)  at  any  givoi  instant  in  time  as  defmed  by  the  user 
from  the  shear  stress  q  and  the  mean  normal  stress  p  using; 

mobilized  4>  =  sin‘‘  (■^)  (57) 

P 

Fw  samples  with  preferred  horizontal  bedding  (0  =  0°),  mobilized  <|>  achieves  higher  peak  values  as  particle 
aspect  ratio  increases.  In  addition,  these  higher  values  are  achieved  at  hi^r  failure  strains.  It  should  be  noted 
that  the  axial  strain  is  in  the  direction  of  Oj  (or  vertical).  For  sarujles  with  60°  bedding,  mobilized  <|>  is  similar 
regardless  of  an  increase  in  particle  aspect  ratio.  No  distinct  peak  is  achieved  in  the  samples  as  their  peak 
strength  and  their  residual  strength  appear  to  be  about  the  same  value. 

Figures  108  and  109  illustrate  the  effect  of  particle  orientation  on  the  progression  of  mobilized  <j>  during  the 
course  of  def(xmati(Ni  for  samples  with  a^  =  1.4  and  3.0,  respectively.  Higher  peak  values  for  mobilized  <|>  are 
achieved  by  assemblages  with  preferred  horizontal  bedding  (0  =  0°).  As  the  bedding  plane  becomes  aligned 
toward  45°  +  <t»/2,  <!>„,„  decreases.  Fot  samples  consisting  of  "Hat"  particles  (a^  =  3.0),  the  difference  in  the 
evoludm  of  mobilized  <t>  due  to  the  orientation  of  the  bedding  plane  is  more  strikmg  as  samples  with  horizontal 
bedding  yield  significantly  higher  values  for  mobilized  <j>  than  those  with  30, 60  or  90°  bedding. 

Similar  trends  in  the  strength  of  the  numerical  san^)les  due  to  particle  aspect  ratio  and  bedding  are  shown 
by  the  evolutiai  of  the  shear  stress  q  and  stress  ratio  C1/02  during  the  course  of  deformation.  Figures  110  through 
113  present  the  results  of  the  biaxial  ccxnpression  test  simulations  for  q,  while  the  results  fw  O1/02  are  presented 
in  Figures  1 14  through  1 17. 
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0  2  4  6  8  10  12 

Axial  Strain  (%) 


♦  a/b=1  -®-a/b=1.1  -*-a/b=1.4  -^a/b=1.7 

^a/b  =  2  ^a/b  =  2.5  ^a/b  =  3 _ 

Figure  106.  Effect  of  particle  aspect  ratio  on  the  mobilized  internal  friction  angle  ({>  for  DEM  samples  with 
preferred  horizontal  bedding  (0  =  0°). 
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Axial  Strain  (%) 


—  a/b=1  -B-a/b=1.1  -*-a/b=1.4  -»-a/b=1.7 

^a/b  =  2  ^a/b  =  2.5  ->-a/b  =  3 _ 

Figure  107.  Effect  of  particle  aspect  ratio  on  the  internal  friction  angle  ^  for  DEM  samples  with  preferred 

60®  bedding. 
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■m-qo  -^30°  —500  -^90«»  I 

Figure  108.  Effect  of  preferred  particle  bedding  0  on  the  mobilized  internal  fricticMi  angle  <|)  for  DEM 
samples  with  a^  =  1.4. 


-m-Qo  30°  -^60°  -^90°  I 

Figure  109.  Effect  of  preferred  particle  bedding  0  on  the  mobilized  internal  friction  angle  <|>  for  DEM 
samples  with  a^  =  3.0, 
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Figure  110.  Effect  of  particle  aspect  ratio  on  the  maximum  shear  stress  q  for  DEM  samples  with  preferred 

horizontal  bedding  (0  =  0°). 


Figure  111.  Effect  of  particle  aspect  ratio  on  the  maximum  shear  stress  q  for  DEM  samples  with  preferred 
60°  bedding. 
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Axial  Strain  (%) 


Figure  112.  Effect  of  piefmed  particle  bedding  6  on  the  maximum  shear  stress  q  for  DEM  samples  with  a/b 

=  1.4. 


Axial  Strain  (%) 


Figure  113.  Effect  of  prefered  particle  bedding  9  cm  the  maximum  shear  stress  q  for  DEM  samples  with  a/b 

=  3.0. 
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—  a/b=1  -B-a/b=1.1  -•-a/b=1.4  -*-a/b  =  1.7 

-^a/b  =  2  -^a/b  =  2.5  —  a/b  =  3 _ 

Figure  114  Effect  of  particle  aspect  ratio  on  the  principal  stress  ratio  a,/<T2  for  DEM  samples  with  preferred 
horizontal  bedding  (0  =  0®). 


Figure  115.  Effect  of  particle  aspect  ratio  cxi  the  principal  stress  ratio  o,/o2  for  DEM  samples  with  preferred 

60®  bedding. 
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a/b=1.4 


-■-0°  -A- 30°  -^60°  -*-90°  I 

Figure  116.  Effect  of  preferred  particle  bedding  0  on  the  principal  stress  ratio  o,/a2  for  DEM  samples  with 

a/b=1.4. 


-w-Qo  ^30°  -*-60°  -^90°  I 

Figure  117.  Effect  of  preferred  particle  bedding  0  on  the  principal  stress  ratio  Oj/oj  for  DEM  samples  with 

a/b  =  3.0. 


100 


8.U  Effect  of  the  Initial  Void  Ratio  on  Strength 

Because  the  surface  of  the  ellipse-shaped  particles  is  smooth  and  continuously  curved,  it  was  very 
difficult  to  f(Hm  numerical  assemblages  with  large  variations  in  their  initial  void  ratio  e^  >^le  maintaining 
preferred  particle  bedding.  Even  when  subjected  to  low  coifining  stresses,  the  loose  assemblages  were  extremely 
unstable  ^  would  immediately  collq>se  into  a  dense  packing.  Various  techniques  were  used  to  consolidate  the 
numerical  samples  to  obtain  loosa  packing  but  most  of  these  packing  schemes  involved  prestressing  the  systems. 
These  techniques  are  described  in  Secdcxi  6.  The  results  of  the  biaxial  shear  simulations  for  these  prestressed 
numoical  sanq)les  are  presented  later  in  Section  8.4.  Some  variation  in  e^  was  achieved  for  the  hand-packed  0 
sample  series  with  zJb  =  2.0  and  0°  bedding.  These  samples  were  not  prestressed  prior  to  the  isotropic 
compression  simulation. 

The  effect  of  the  initial  void  ratio  on  for  DEM  samples  consisting  of  particles  with  a/b  =  2.0  and 
preferred  0°  bedding  is  shown  in  Figure  118.  Samples  with  denser  packing  exhibit  higher  strength  =  47®) 
than  those  with  looser  packing  (<{>„,„  =  40®).  The  results  for  the  numerical  assemblages  are  similar  to  the  those 
reported  by  Holubec  and  D'Appolonia  (1973)  for  different  sands  and  glass  beads  shown  in  Figure  1 19. 

The  effect  of  the  initial  void  ratio  e„  is  similar  on  the  peak  shear  stress  q  and  the  stress  ratio  q  /q  as 
shown  in  Figures  120  and  121.  Higher  strength  is  exhibited  by  the  denser  sample  (q„,„  =  2.7  kPa,  (ofi/o2)iMK  “ 
6.3)  than  the  looser  sample  (q^  =1.8  kPa,  (0^02)0^  =  4.7).  The  evolution  of  mobilized  <|),  q  and  (t,/c2  for  the 
numerical  samples  with  variations  in  e^  are  shown  in  Figures  122  through  124,  respectively. 


dense  Initial  Void  Ratio  loose 


Figure  118  Effect  of  die  initial  void  ratio  on  the  peak  friction  angle  for  DEM  samples  with  a/b  =  2.0  and 
0  =  0®. 
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Figure  119.  Effect  of  the  relative  density  on  the  internal  friction  angle  (j>  for  different  granular  materials 
(from  Holubec  and  d’Appolonia  1973). 


Initial  Void  Ratio  loose 


Figure  120.  Effect  of  the  initial  void  ratio  on  the  peak  shear  stress  for  DEM  samples  with  a^  =  2.0 

and  0  =  0®. 
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Figure  121.  Effect  of  the  initial  void  ratio  cm  the  peak  principal  stress  ratio  (oi/cT2)m<«  for  DEM  samples  with 

a/b  =  2.0  and  0  =  0°. 
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♦  6  =  0.159 


e  =  0.166  -*-6  =  0.176 


Figure  122.  Effect  of  the  initial  void  ratio  e  on  the  mobilized  internal  friction  angle  <|)  for  DEM  samples  with 

a/b  =  2.0  and  0  =  0®. 
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a^  =  2.0 


-—e  =  0.159  -^e  =  0.166  —e  =  0.176  | 

Figure  123.  Effect  of  the  initial  void  ratio  e  rai  the  maximum  shear  stress  q  Cw  DEM  samples  with  a/b  =  2.0 

and  0  =  0®. 


-■-6  =  0.159  -^e  =  0.166  —e  =  0.176  | 

Figure  124.  Effect  of  the  initial  void  ratio  e  on  the  principal  stress  ratio  cjci  for  samples  with  a^  = 
2.0  and  0  =  0° 
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8.1 .4  Effect  of  a  Random  Orientation  Fabric 


The  orientation  fabric  of  the  DEM  tests  presented  thus  far  have  been  extremely  anisotropic.  Particles 
werenumoically  generated  with  preferred  orientaticms  of 0, 30, 60  and  90°  within  a  standard  deviation  of  15°. 
Additional  sanples  were  generated  with  a  random  distribution  of  particle  orientations  creating  a  nearly  isotropic 
(xioitatiai  fabric.  The  samples  with  random  particle  bedding  were  significantly  looser  than  those  with  aligned 
particles.  For  example,  the  respective  initial  void  ratios  for  Oval  2  (a/b  =  1.4)  and  Oval  7  (a/b  =  3.0)  samples 
with  randcnn  bedding  were  0.204  and  0.236,  while  those  with  preferred  horizontal  bedding  (0  =  0°)  averaged 
0. 174  and  0. 170.  Presumably,  the  nature  of  randcxn  bedding  produces  much  looser  packing  as  more  side  to  edge 
ccmtacts  exist  betweoi  the  randnnly  oiented  particles.  Samples  with  anisotropic  fabrics,  on  the  other  hand,  tend 
to  be  denser  as  more  side  to  side  contacts  exist  between  the  aligned  particles. 

Figures  125  to  127  show  the  effect  of  particle  aspect  ratio  on  the  peak  internal  friction  angle  (t>^,  peak 
shear  stress  and  peak  stress  ratio  (?  /?  for  samples  with  preferred  horizontal  bedding  and  random 
bedding.  While  the  overall  strength  of  the  assemblages  with  both  preferred  and  random  bedding  increases  as 
particle  aspect  ratio  increases,  the  peak  values  are  significantly  lower  for  samples  with  random  bedding.  Oval 
2  and  Oval  7  samples  with  random  bedding  yield  respective  values  for  <|)„^  of  26  and  3 1°,  while  those  with 
prefoied  horizontal  bedding  (0=  0°)  yield  average  values  of  36  and  55°.  This  is  due  in  part  to  the  looser  packing 
in  the  random  samples. 


■  0°  bedding  -a-  random  bedding  \ 

Figure  1 25.  Effect  of  particle  aspect  ratio  on  the  peak  friction  angle  for  DEM  samples  with  preferred 

horizontal  bedding  (6  =  0°)  and  random  bedding. 
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■  0°  bedding _ -e-  random  bedding  \ 

Figure  126,  Effect  of  particle  aspect  ratio  on  the  peak  shear  stress  for  DEM  samples  with  preferred 

horizontal  bedding  (6  =  0®)  and  random  bedding. 


■  0°  bedding _ -b-  random  bedding  { 

Figure  127.  Effect  ofparticle  aspect  ratio  cm  the  peak  principal  stress  ratio  (oj/aj)^  for  DEM  samples  with 

preferred  horizontal  (0  =  0®)  and  random  bedding. 
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8.2  Influence  of  Particle  Shape  on  Volumetric  Behavior 

The  volumetric  befaavio'  of  granular  matoials  is  also  strongly  influenced  by  particle  shape.  Angular  particles 
interlock  m(xe  which  facilitates  more  compression  during  the  initial  stages  of  shear.  The  tendency  for 
assemblages  with  angular  particles  to  contract  more  initially  was  observed  in  the  volumetric  strain  behavior  of 
both  Kraishi  et  al's  e?q)aimental  tests  and  the  ellipse-based  DEM  simulations  presented  in  Section  7.  However, 
as  defamation  continues,  the  angular  particles  ultimately  dilate  more  than  the  more  roimded  particles  when  the 
orientation  of  the  preferred  bedding  plane  is  nomal  to  the  direction  of  the  major  principal  stress.  The  data 
caiq>arison  also  showed  that  whoi  the  bedding  plane  is  aligned  toward  the  predicted  Mohr-Coulomb  failure  plane 
(45° + (ti/2),  the  angular  particles  dilate  less  flian  the  more  rounded  particles  during  biaxial  shear.  This  difference 
in  the  dilatant  behavior  due  to  particle  aspect  ratio  and  bedding  is  examined  further  in  Section  9. 

8.2.1  Effect  of  Particle  Aspect  Ratio 

The  effect  of  particle  aspect  ratio  on  the  volumetric  strain  behavior  for  DEM  samples  with  preferred 
hoTzootal  bedding  (0  =  0°)  is  shown  in  Figure  128.  As  shown  previously,  the  sign  conventions  for  the  volumetric 
strain  plots  shown  are  opposite  fixm  those  nomalty  associated  with  contraction  and  dilation.  Though  the  vertical 
axial  strains  shown  are  positive  ftx'  compression,  Ae  volumetric  strains  shown  are  negative  for  compression  and 
positive  for  dilation. 

As  particle  aspect  ratio  increases,  the  assemblages  experience  longer  periods  of  initial  contraction  but 
ultimately  achieve  higher  dilatant  strains  over  the  course  of  defamation.  At  larger  axial  strains,  the  samples 
achieve  a  peak  after  which  a  state  of  constant  volume  strain  exists.  This  coitractant  and  dilatant  behavior  is 
typical  for  dense  soils  and  other  granular  materials. 
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Figure  1 28.  Effect  of  particle  aspect  ratio  on  the  volumetric  strain  for  DEM  samples  with  0  =  0° 
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Similar  behavior  is  observed  in  the  evolution  of  the  void  ratio  e  during  the  course  of  stress*induced 
(tftfnrmatinn  as  illustrated  in  Figure  129.  Priw  to  shear,  the  angular  particles  form  denser  assemblages  than  the 
more  rounded  particles  as  they  are  able  to  interlock  more.  Dunng  shear,  the  samples  with  angular  particles  and 
prefarcd  horizoital  bedding  dilate  more.  The  relative  increases  in  e  from  its  initial  value  is  significantly  higher 
for  samples  with  angular  particles.  For  example,  a  15%  increase  in  the  voids  is  yielded  by  samples  with  a/b  — 
1.0  (fiom  eo = 0.20  to  e  =  0.23)  and  56%  in  those  with  a^  =  2.5  (from  0. 16  to  0.25).  After  approximately  10% 
axial  strain,  the  ass^blages  with  rounded  particles  approach  a  state  of  constant  volume  strain,  while  those  with 
more  angular  particles  continue  to  dilate. 


♦  a/b=  1 

-e- a/b  =1.1 

a/b  =  1.4  a/b  =  1.7 

II 
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a/b  =  2.5 

-•-a/b  =  3 

Figure  129.  Effect  of  particle  aspect  ratio  on  the  void  ratio  for  DEM  samples  with  preferred  horizontal 
bedding  during  biaxial  shear. 


8.2.2  Effect  of  Particle  Aspect  Ratio  and  Orientation 

The  effect  of  particle  aspect  ratio  on  the  overall  deformation  behavior  is  significantly  different  for  samples 
with  preferred  bedding  planes  of  30, 60  or  90®.  As  particle  aspect  ratio  increases,  the  assemblages  contract 
IrmgfT  and  Hilate  less  than  the  san^les  with  horizmtal  bedding  (0  =  0°).  The  volumetric  strain  for  DEM  samples 
with  60®  bedding  is  shown  in  Figure  130.  Though  the  samples  have  similar  initial  densities  (or  void  ratios),  the 
volumetric  strain  behavior  exhibited  is  similar  to  that  typically  exhibited  by  loose  materials.  As  the  bedding  plane 
becomes  aligned  toward  the  orientation  of  the  Mohr-Coulomb  failure  plane  (45°  +  <t>/2),  dilation  actually 
decreases  with  increased  particle  aspect  ratio.  This  trend  is  due  to  the  coaxial  alignment  of  the  interlocking  plane 
(ot  bedding  plane)  with  45®  +  <|)/2  [8].  When  these  planes  bewane  coaxial,  assemblages  with  higher  degrees  of 
intffrlfvrlfing  deform  more  as  a  "rigid  body"  thereby  exhibiting  significantly  less  overall  volumetric  strain.  This 
behavior  is  explained  further  in  Chapter  9. 
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Figure  1 30.  Effect  of  particle  aspect  ratio  on  the  volumetric  strain  for  DEM  samples  with  6  =  60® 

Figure  131  shows  the  effect  of  prefan*ed  particle  orientation  on  the  volumetric  strain  behaviOT  for  samples 
consisting  of  "flat"  particles  (a/b  =  3.0).  Significantly  more  dilation  occurs  when  the  bedding  plane  is  aligned 
iKxmal  to  die  majcx'  principal  stress  (0  =  0®)  than  when  the  bedding  plane  (or  interlocking  plane)  is  30, 60  or  90®. 


Figure  131.  Effect  of  particle  bedding  on  the  volumetric  strain  for  DEM  samples  with  a/b  =  3.0. 
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While  dififooices  in  the  deformation  behavior  of  granular  materials  is  normally  associated  with  differences 
in  the  sample's  relative  density,  the  numoical  assemblages  possess  approximately  the  same  initial  void  ratios  (e^ 
=  0. 174, 0. 186, 0. 178  and  0. 178  fw  0, 30, 60  and  90°  bedding,  respectively).  Similar  trends  are  observed  in  the 
evolution  of  the  void  ratio  shown  in  Figures  132  and  133. 


II 

-B-a/b  =  1.1 

-«h-a/b  =  1.4  -A-a/b  =  1.7 
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a/b  =  2.5 
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Figure  132.  Effect  of  particle  aspect  ratio  on  the  void  ratio  for  DEM  samples  with  preferred  60°  bedding 
during  biaxial  shear. 


-»-0°  -^30°  -^60°  -^90°  I 

Figure  133.  Effect  of  particle  orientation  on  the  void  ratio  for  DEM  samples  with  a^  =  3.0. 
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The  angle  of  dilatancy  5  rq)resents  the  overall  volumetric  strain  rate  for  the  DEM  assemblages  from  the  start 
of  shear  (e,  0%)  to  the  strain  at  which  AV\y,  is  maximum.  The  following  equation  is  used  to  compute  the 

dilatancy  angle  for  the  rectangular  DEM  samples  under  biaxial  conq)ressi(»i; 


sin6  = 


dz^-dz^ 


(58) 


whoe  de„  dy^^,  de,  and  dcj  are  the  changes  in  volumetric  strain,  shear  strain,  axial  (vertical)  strain  and  lateral 
(hcxizmtal)  strain,  respective^.  The  strain  incremoits  are  positive  for  crnnpression  and  negative  for  dilation  and 
are  taken  over  the  same  time  increment. 


The  effects  of  particle  aspect  ratio  and  the  orientation  of  the  bedding  plane  on  5  for  the  DEM  samples 
are  shown  in  Figures  134  and  135.  Assemblages  consisting  of  round  particles  (a/b  =  1.0)  yield  values  for  5  of 
7.5°,  while  samples  with  angular  particles  and  horizontal  bedding  (0  =  0°)  yield  values  for  6  of  approximately 
15°.  While  5  increases  with  increased  aspect  ratio  for  samples  with  6  =  0°,  it  decreases  significantly  with 
increased  aspect  ratio  fix'  sanq)les  with  0  =  30, 60  and  90°.  The  difference  in  5  due  to  particle  orientation  is  small 
fix  sanq)les  coisisting  of  rounded  particles,  but  as  aspect  ratio  increases  and  the  orientation  fabric  becomes  more 
anisotropic,  the  difference  is  much  greater.  The  angle  of  dilatancy  is  lowest  for  samples  consisting  of  angular 
particles  with  60°  bedding  as  they  approach  5  =  0°  indicating  almost  no  dilation  behavior. 

The  defomaticxi  of  the  particulate  assemblages  with  0  =  0°  can  be  related  to  the  strength  as  shown  in  Figure 
136.  Hitler  strength  is  achieved  by  san:q>les  that  experience  more  overall  dilation  during  the  deforming  process. 


■  0“  A  30°  •  60°  ▼  90° 


Figure  134.  Effect  of  particle  aspect  ratio  on  the  dilatancy  angle  for  DEM  samples  with  preferred  particle 
bedding. 
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Figure  135.  Effect  of  particle  bedding  on  the  dilatancy  angle  for  DEM  samples  with  variations  in  particle 
aspect  ratio. 


Figure  136.  Effect  of  the  dilatancy  angle  cm  the  peak  internal  friction  angle  for  DEM  samples  with  a/b  —  2.0 
and  0  =  0°. 
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8.2  J  Effect  of  Void  Ratio 


The  effect  of  the  initial  void  ratio  e^  on  the  volumetric  strain  of  DEM  samples  with  a/b  =  2.0  and 
preferred  horizcmtal  bedding  (0  =  0®)  is  shown  in  Figure  137.  The  sample  with  looser  packing  (e,,  =  0.176) 
exhibits  less  dilatant  behavior  than  the  denser  sample  (e^  =  0. 158).  Similarly,  Figure  138  shows  the  evolution 
of  the  void  ratio  during  the  course  of  deformation  for  the  same  numerical  samples.  While  each  sample  exhibits 
"dense"  bdiavior,  the  doiser  sanq)le  (e^ = 0.158)  dilates  more  than  the  looser  sample.  This  is  typical  defoimatirai 
behavior  fcH*  "real"  granular  materials  with  variations  in  their  initial  relative  density. 

The  effect  of  the  void  ratio  on  the  dilatancy  angle  5  for  DEM  samples  with  a/b  =  2.0  and  preferred 
hcxizmtal  bedding  (6  =  0°)  is  shown  in  Figure  139.  Again,  the  tendency  for  the  denser  sample  to  dilate  more  than 
the  less  dense  sample  is  apparent.  The  assemblage  with  e^  =  0. 158  yields  a  value  for  5  of  15°,  while  for  ^  = 
0. 176, 6  is  approximately  7  or  8°. 
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Figure  137.  Effect  of  void  ratio  on  the  volumetric  strain  for  DEM  samples  with  a/b  =  2.0  and  0  =  0°. 
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Figure  138.  Evoluticm  of  die  void  ratio  e  for  DEM  san^les  with  a/b  =  2.0, 0  =  0°  and  variations  in  the  initial 
void  ratio. 


Initial  Void  Ratio  loose 


Figure  139.  Effect  of  the  initial  void  ratio  on  the  dilatancy  angle  for  DEM  samples  with  a/b  =  2.0  and  0  = 
0°. 
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8.2.4  Effect  of  an  Isotropic  Orientation  Fabric 

The  progressions  of  the  volumetric  strain  for  numerical  assemblages  with  random  bedding  are  shown  in 
Figure  140..  The  randcxn  sanqrle  assisting  of  round  particles  is  approximately  isotropic  and  experiences  almost 
no  volume  change  during  the  course  of  defcxmaticat  As  particle  aspect  ratio  increases,  however,  the  assemblages 
contract  and  dilate  slightly  more.  Because  the  fabrics  are  looser  and  significantly  less  anisotropic  than  the 
assemblages  with  preferred  bedding,  the  random  samples  contract  and  dilate  less.  The  small  increase  in  the 
dilatant  bdiavkxr  fix  the  random  assemblages  consisting  of  flatter  particles  is  due  to  an  increase  in  interlocking. 

Figure  141  illustrates  the  effect  of  particle  aspect  ratio  on  the  angle  of  dilatani^  5  for  DEM  samples  with 
preferred  horizontal  and  randcxn  bedding.  The  samples  with  more  isotropic  fabrics  (random  bedding)  exhibit 
significantly  less  dilatant  bdiavior  than  the  samples  with  preferred  bedding.  Though  the  numerical  assemblages 
with  preferred  bedding  yield  values  for  5  that  increase  with  particle  aspect  ratio,  the  samples  with  random  bedding 
remain  at  approximately  the  same  value. 

The  fabric  of  assonblages  consisting  of  perfectly  round  particles  (a^  =  1.0)  are  more  isotropic  since  round 
particles  have  no  orientation.  Consequently,  the  difference  in  5  for  these  samples  is  primarily  attributed  to  the 
difference  in  their  initial  void  ratios. 


•*-a/b=1.0  ^a/b  =  2.0  -^a/b  =  3.0  | 

Figure  1 40.  Effect  of  particle  asp)ect  ratio  on  the  volumetric  strain  for  DEM  samples  with  random  bedding. 
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Particle  Angularity  (a/b) 


■  0^  bedding _ random  bedding  | 

Figure  141 .  EflFect  of  particle  aspect  ratio  on  the  dilatancy  angle  for  DEM  samples  with  preferred  horiztmtal 

bedding  (0  =  0°)  and  randcmi  bedding. 


8.2.5  Deformation  Mechanisms 

Previous  micranechanical  investigations  have  shown  that  particle  rolling  may  be  the  dominant  deformation 
mechanism  and  primary  cause  for  failure  of  particulate  systems  (Oda  et  al  1983,  Ting  et  al  1993).  Numerical 
models  using  circles  or  sphaes  for  particle  sh^)es  have  ctaisistently  yielded  lower  overall  strength  compared  with 
real  sands  because  die  rcxmd  particles  toid  to  roll  excessively  during  the  assemblage's  deformation.  For  samples 
with  round  particles  (a/b  =  1.0)  in  the  current  study,  the  peak  angle  of  internal  friction  ranges  from  a  low 

of  19°  in  looser  assemblages  to  a  high  of  26°  in  those  with  denser  packing.  Since  <j>TO«x  o**  average  is  lower  than 
the  interparticle  sliding  friction  angle  <|>^  of  26°,  it  is  likely  that  the  low  strength  exhibited  by  assemblages  with 
roimd  particles  is  due  to  particle  rolling.  The  primaiy  advantage  to  using  angular  shapes  such  as  the  ellipse  is 
that  particle  interlocking  can  occur  which  significantly  reduces  the  negative  effects  of  particle  rolling  and  more 
accurately  simulates  the  micrcmiechanical  behavior  of  real  particulate  systems. 

The  ellipse-based  DEM  model  tracks  the  motion  of  each  particle  during  the  course  of  the  simulation  and 
decanposes  the  relative  contact  velocities  between  colliding  particles  into  translation  and  rotation  ounponents. 
The  total  rotation  and  translation  (xunponents  of  the  relative  contact  velocities  are  further  decomposed  into 
normal  and  tangential  components.  These  incremental  quantities  are  summed  over  every  contact  within  the 
assemblage  to  assess  the  overall  relative  importance  of  particle  motion  due  to  rolling  and  sliding  during  the 
assemblage's  deformation.  Very  little  net  contact  motion  occurs  in  the  subdomains  within  the  assemblages  where 
rigid  bo(fy  motitm  exists.  Consequently,  net  contact  deformations  are  used  instead  of  total  contact  deformations 
to  effectively  exclude  the  rigid  body  motion  from  the  contact  deformation  statistics.  Details  of  how  the  model 
decomposes  these  contact  deformations  have  been  presented  elsewhere  (Ting  et  al  1993). 
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Figure  142  shows  the  effect  of  particle  aspect  ratio  on  the  relative  importance  of  rotation  and  translation  cm 
the  total,  normal  and  tangential  particle  motion  for  the  DEM  samples  with  preferred  bedding  of  0°.  For  the 
numerical  samples  consisting  of  round  particles  (a/b  =1.0),  the  net  total  particle  rotation  is  over  2  times  higher 
than  particle  translation.  No  net  norm^  rotation  occurs  for  the  round  particles  since  their  rotation  caimot  exert 
eccentric  tKxmal  faces  on  nei^boring  particles.  Instead,  contact  motion  due  to  rotation  between  round  particles 
is  entirely  tangential.  As  particle  aspect  ratio  increases,  however,  the  relative  importance  of  particle  motion  due 
to  rotatioi  and  translation  converges  toward  unity.  This  is  true  for  the  total,  normal  and  tangential  components 
of  particle  motion.  This  suggests  that  angular  particles  exert  eccentric  normal  forces  on  neighboring  particles 
in  contact  thereby  significantly  reducing  the  negative  effects  of  particle  rolling. 

Similar  results  are  observed  in  Figure  143  for  DEM  samples  with  preferred  30°  bedding.  The  similarities 
indicate  that  the  relative  importance  of  particle  motion  due  to  rotation  and  translation  are  independent  of  the 
bedding  plane's  orientation  with  respect  to  the  orientation  of  the  principal  stress  a,. 


Figure  142.  Effect  of  particle  aspect  ratio  a/b  on  the  relative  importance  of  deformation  mechanisms  fix' 
DEM  samples  with  preferred  0°  bedding. 
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Figure  143.  Effect  of  particle  aspect  ratio  a/b  on  the  relative  importance  of  deformation  mechanisms  for 

DEM  samples  with  preferred  30°  bedding. 


83  Influence  of  Particle  Shape  on  the  Coordination  Number 

The  rnnrHinarifin  number  €„  is  the  average  numbCT  of  ccmtacts  per  particle  and  represents  the  local  variability 
of  the  assemblage's  pwosity  (w  void  ratio).  Studies  have  shown  that  samples  with  higher  coordination  numbers 
told  to  exhibit  higher  strength  as  more  particle  contacts  produce  a  denser  and  more  stable  sample  (Feda  1982). 
The  effect  of  particle  aspect  ratio  cm  the  initial  cocmdination  number  for  the  numerical  samples  with  0, 30, 60  and 
90°  bedding  is  shown  in  Figure  144.  The  initial  C„  increases  with  increasing  particle  aspect  ratio  and  is 
apparently  unaffected  by  the  orientaticm  of  the  overall  preferred  bedding  plane. 

During  the  course  of  stress-induced  deformation,  however,  the  orientation  of  the  bedding  plane  as  well  as 
particle  aspect  ratio  have  a  significant  effect  on  the  evolution  of  the  coordination  number.  Figures  145  and  146 
illustrate  the  effect  of  particle  aspect  ratio  on  C„  for  samples  with  0  and  60°  bedding,  respectively. 
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Coordination  Number 


Particle  Angularity  (a/b) 


Figure  144.  Effect  of  particle  aspect  ratio  on  the  initial  coordination  number  for  DEM  samples  with 
preferred  bedding  6. 


Figure  145.  Effect  of  particle  aspect  ratio  on  the  coordination  number  for  DEM  samples  with  0  =  0°. 


As  the  samples  with  horizontal  bedding  (0  =  0®)  are  sheared,  relatively  little  change  in  C„  is  observed  for 
those  samples  with  a/b  of  1.4  or  less.  For  samples  consisting  of  more  angular  particles  (a/b  >  1.4)  and  0° 
bedding,  however,  C„  decreases  from  its  high  initi^  value  to  a  significantly  lower  value  that  is  similar  to  that  of 
the  samples  coisisting  of  more  rounded  particles.  For  example,  C„  for  samples  c(Hisisting  of  "flat"  particles  (a/b 
=  3.0)  decreases  from  a  value  of  3.8  prior  to  shear  to  a  value  of  approximately  3.1  at  10%  axial  strain. 
Collectively,  the  san:q)les  with  prefared  horizontal  bedding  ^pear  to  be  converging  to  a  similar  residual  C„  value 
of  ^proximately  2.8  as  the  particles  beccxne  fiiUy  mobilized.  In  another  DEM  study  in  which  randmnly  oriented 
and  i^orm  siz^  ellipse-sluqied  particles  woe  i^,  Rothenburg  and  Bathurst  (1991)  found  that  C„  approached 
a  "critical"  value  of  approximately  3.5.  The  fabric  of  their  numerical  samples,  however,  were  mwe  isotropic. 

Cn  for  the  samples  with  60®  bedding,  however,  does  not  tend  to  converge.  Instead,  for  those  samples 
cmisisting  of  more  angular  particles,  the  coordination  number  quickly  increases  from  its  initial  value  and  then 
remains  relatively  constant  during  the  course  of  deformation.  For  example,  C„  for  samples  consisting  of  flatter 
particles  (a/b  =  3.0)  inaeases  frem  its  initial  value  of  3.8  to  a  value  of  4.2  at  approximately  4%  axial  strain  and 
thereafter  remains  about  the  same. 

The  effect  of  particle  bedding  on  the  evolution  of  the  coordination  number  for  DEM  samples  with  a/b  =  1. 1 
and  3.0  is  shown  in  Figures  147  and  148,  respectively.  For  samples  consisting  of  rounded  particles  (a/b  =1.1), 
the  effect  of  the  bedding  plane  is  rather  insignificant  as  similar  values  for  C„  are  obtained  for  the  numerical 
samples  with  variations  in  particle  bedding.  During  biaxial  shear,  C„  remains  approximately  the  same  for 
samples  with  a/b  =  1.1.  The  evolution  of  C„  is  strongly  influenced,  however,  as  the  particles  become  more 
angular.  For  samples  consisting  of  flatter  particles  (a/b  =  3.0)  with  preferred  horizontal  bedding  (0  =  0®),  the 
coordination  number  decreases  significantly  during  the  course  of  deformation.  The  decrease  in  C„  is  less  for 
sanples  with  30°  bedding,  while  for  samples  with  60  and  90°  bedding,  C„  actually  increases.  Similar  decreases 
in  the  cocxdination  number  have  been  reported  by  Rothenburg  and  Bathurst  (1991,  Chen  and  Ishibashi  (1990). 


—  a/b=  1 

-*-a/b=  1.1 

^  a/b  =  1.4  a/b  =  1.7 

CM 

II 

a/b  =  2.5 

-^a/b  =  3 

Figure  146.  Effect  of  particle  aspect  ratio  on  the  coordination  number  for  DEM  samples  with  0  -  60®. 
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Effect  of  particle  bedding  on  the  coordination  number  for  DEM  samples  with  a/b  =  1. 1. 


Figure  148.  Effect  of  particle  bedding  on  the  coordination  number  for  DEM  samples  with  a^  =  3.0. 


This  is  in  agreement  with  the  deformation  behavior  of  the  particulate  assemblages  presented  previously. 
When  subjected  to  biaxial  shear,  sanq>les  coisisting  of  angular  particles  with  preferred  horizontal  bedding  dilate 
more  than  those  with  30, 60  or  90®  bedding.  As  a  result,  more  contacts  are  broken  between  particles  as  the 
assonblage  dilates,  while  mwe  ccmtacts  are  made  between  particles  with  60  and  90°  bedding  as  the  assemblage 
contracts.  This  trend  is  not  as  distinct  in  samples  consisting  of  rounded  particles,  however,  because  less 
differences  exist  in  the  dilatant  behavior  due  to  particle  bedding. 

The  effect  of  particle  aspect  ratio  cm  the  cocxdinaticm  numbo'  fix  numerical  samples  with  preferred  horizontal 
hfAting  (0  =  0°)  and  random  bedding  is  shown  in  Figure  149.  Figure  144  showed  no  significant  differences  in 
C„  for  samples  with  preferred  particle  bedding  of  0, 30, 60  and  90°.  Similarly,  no  significant  differences  in  C„ 
are  evident  fcx  assemblages  with  preferred  and  random  bedding.  This  is  unexpected  since  the  relative  densities 
of  the  assemblages  are  quite  different.  Because  C„  is  generally  regarded  as  being  closely  associated  with  the 
relative  density,  it  is  expected  that  the  coordination  number  for  the  looser  samples  with  random  bedding  to  be 
lower  than  few  those  with  preferred  bedding.  However,  Chen  and  Ishibashi  (1990)  suggest  in  their  numerical 
simulations  using  spherical  particles  that  isotropic  fabrics  have  the  highest  ccxirdination  number. 


■  0°  bedding _ -a-  random  bedding  | 

Figure  149.  Effect  of  particle  aspect  ratio  on  the  initial  coordination  number  for  DEM  samples  with 
preferred  horizontal  bedding  (0  =  0°)  and  random  bedding. 
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8.4  Effect  of  Sample  Formation 

The  majority  of  the  numerical  assemblages  were  formed  by  individually  placing  the  particles  in  their  initial 
locatkxis  forming  densely  padced  sanq)les  (0  sample  saies).  In  an  effort  to  obtain  samples  with  more  variations 
in  the  initial  void  ratio  e„  viiile  maintaining  control  of  particle  bedding,  other  packing  techniques  were  used  to 
form  the  L,  FD  and  FL  sample  series.  A  wider  range  for  e^  was  achieved  using  the  different  packing  schemes, 
however  the  ranges  for  individual  packing  methods  were  still  rather  narrow.  For  example,  the  L  sample  series 
yielded  values  for  Oo  ranging  from  0. 185  to  0.21 1,  while  for  the  FD  series,  ranged  from  0. 139  to  0. 149. 

The  effect  of  the  initial  void  ratio  and  the  different  packing  methods  on  the  peak  internal  friction  angle 
for  sanq)les  with  a/b  =  2.0  and  6  =  0°  is  shown  in  Figure  150.  Though  e^  varies  significantly. for  the  combinaticm 
of  sample  series,  the  magnitude  of  <!>„„  does  not  vary  accordingly.  Instead,  thq,ej4  relationship  for  the 
individual  sample  series  shifts  aloig  the  e^  axis  while  ({^nax  i^ains  relatively  unchanged.  These  results  illustrate 
the  extrone  sensitivity  of  the  ellipse-based  DEM  model  to  the  stress  history  of  the  samples  during  their  formation. 

The  progression  of  the  maximum  shear  stress  q  and  volumetric  strain  for  individual  numerical 
samples  with  a/b  =  2.0  and  0°  bedding  is  shown  in  Figures  151  and  152.  Assemblages  that  were  prestressed 
during  their  consolidation  (FD51 14  and  L51 12)  generally  exhibit  stiffer  responses  to  shear,  while  the  hand- 
padctd  sample  (05 1 14)  >riudi  was  not  prestressed  exhibits  a  softer  response.  Sample  FL5 1 14  exhibits  a  softer 
response  similar  to  the  hand-packed  sample  suggesting  that  its  initial  stress  fabric  prior  to  shear  is  less 
anisotix^ic  than  samples  FD5 114  and  L5 1 1 2.  Presumably,  the  more  anisotropic  stress  fabric  produces  a  more 
stable  sanq)le  initially  and  Inittle  stress-strain  behavior,  while  less  stress  fabric  anisotropy  produces  a  less  stable 
sanq>le  and  plastic  b^vicx.  Similar  stiff  responses  to  shear  due  to  prestressing  was  observed  by  Ishibashi  and 
Agarwal  (1991)  in  their  numerical  simulations  using  spheres. 


Initial  Void  Ratio  loose 


■  O  series  •  L  series  a  fD  series  a  fL  series 


Figure  150.  Effect  of  the  initial  void  ratio  <mi  the  peak  friction  angle  for  DEM  samples  (a/b  =  2.0  and 

0  =  0®)  formed  differently 
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-■-05114  -iir-FD5114  -^FL5114  —LS1 12  | 

Figure  151.  Efifect  of  the  initial  void  ratio  cm  the  shear  stress  for  DEM  samples  (a/b  =  2.0  and  9  =  0®) 

formed  differently. 


-■-05114  -■-FDSIU  -■-FL5114  — -L5112  | 

Figure  152.  Effect  of  the  initial  void  ratio  on  the  volumetric  strain  for  DEM  samples  (a/b  =  2.0  and  9  = 

0®)  formed  differently. 
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8.5  Effect  of  Confining  Stress  Oj 

The  data  presetited  thus  far  was  obtained  from  samples  that  were  isotropically  compressed  to  a  mean  stress 
p  of  0.7  kPa  and  then  biaxially  sheared  while  laterally  iq)plying  a  constant  (xmfming  stress  Oj  of  0.7  kPa. 
Additioial  san:q)les  with  particle  aspect  ratios  a/b  oil  A  and  2.0  and  with  preferred  horizontal  bedding  (6  =  0°) 
were  isotropically  ccmipressed  to  a  p  of  1.4  and  2.1  kPa  and  then  sheared  with  =  1.4  and  2.1  kPa  applied 
laterally.  This  Section  presents  the  strength  and  deformation  results  of  these  additional  test  simulations  with 
variations  in  Oj. 

Figure  153  shows  the  effect  of  variaticxis  in  02  on  the  shear  stress  q  for  the  same  numerical  assemblage  with 
a/b  =  1 .4  and  6  =  0°.  An  inoease  in  the  shear  resistance,  hence  the  strength,  is  achieved  by  the  assemblage  with 
increased  Cj-  The  effect  is  similar  for  a  DEM  sample  with  a/b  =  2.0  and  0  =  0®  as  shown  in  Figure  154.  A 
c(xiq>aris(m  of  the  results  shown  in  Figures  153  and  154  shows  the  effect  of  particle  angularity  oa  the  strength 
of  the  assemblages.  For  each  the  san:q>le  consisting  of  angular  particles  (a/b  =  2.0)  yields  significantly  higher 
peak  values  for  q  than  that  with  the  m«e  rounded  particles  (a/b  =  1.4).  For  example,  assemblages  with  a/b  =  1.4 
and  2.0  that  are  confmed  laterally  with  Oj  =  0.7  kPa  yield  values  for  of  1.5  and  2.7  kPa,  respectively,  while 
those  with  O2  =  2. 1  kPa  yield  values  for  of  4.0  and  6.5  kPa. 

The  maximum  shear  stress  q  is  plotted  against  the  mean  stress p  in  Figure  155  for  the  numerical  sample  with 
a/b  =  1 .4  and  in  Figure  156  for  the  sanqile  with  a/b  =  2.0.  The  stress  paths  are  oriented  at  45°  from  the  horizontal 
as  they  are  for  "real"  granular  materials.  The  angle  of  intonal  friction  ({)  is  a  material  property  and  does  not  change 
significantib'  ^  ^  confining  stress  02  is  varied.  The  Mohr-Coulomb  failure  envelope  defines  (|)  as  illustrated  in 
Figures  155  and  156.  The  effect  of  particle  angularity  on  <t>  is  evident  as  the  more  angular  particles  with  a/b  = 
2.0  yield  a  higher  average  value  for  <)»  than  the  more  rounded  particles  with  a/b  =  1.4  (<t>  =  47.0  and  38.4°, 
respectively). 


Figure  152.  Effect  of  the  confining  stress  02  on  the  shear  stress  q  for  a  DEM  sample  with  a/b  =  1.4  and  0° 

bedding. 
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Figure  154.  Effect  of  the  confining  stress  Oj  on  the  shear  stress  q  for  a  DEM  sample  with  a^  =  2.0  and  0° 

bedding. 


Figure  155.  M<Ar  circles  for  a  DEM  sample  (a/b  =  1.4, 0  =  0°)  with  confining  stresses  of  0.7, 1.4  and  2. 1 
kPa  and  the  resulting  Mohr-Coulomb  failure  envelope. 
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Figure  156.  Mdir  circles  for  a  DEM  sample  (a/b  =  2.0, 0  =  0°)  with  confining  stresses  of  0.7, 1.4  and  2. 1 
kPa  and  the  resulting  Mohr-Coulomb  failure  envelope. 


The  effect  of  Cj  oti  the  volumetric  strain  AVA^^  is  shown  in  Figures  157  and  158.  The  DEM  results  show 
that  the  san^le  with  the  hitter  confining  stress  dilates  less  throughout  the  course  of  deformation.  This  behavior 
is  typical  f(^  soil  and  other  granular  materials.  A  comparison  of  Figures  157  and  158  illustrate  the  effect  of 
particle  angularity  c»  AVA^^.  The  assemblage  with  angular  particles  (a^  =  2.0)  dilates  more  than  that  with  more 
rounded  particles. 
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Figure  157. 


Figure  158. 


Effect  of  the  coifining  stress  02  0“  the  volumetric  strain  for  a  DEM  sample  with  a/b  =  1 .4  and 
0®  bedding. 
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Effect  of  the  confining  stress  CTj  on  the  volumetric  strain  for  a  DEM  sample  with  a/b  —  2.0  and 
0®  bedding. 
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9.  DISCUSSION  OF  DEM  RESULTS 


The  results  of  the  biaxial  compression  test  simulations  presorted  in  Section  8  show  that  the  overall 
medranical  behavior  of  granular  materials  is  strongly  influenced  by  the  shape  of  the  particles.  Higher  strength 
and  more  dilatioi  are  exhibited  by  the  numerical  assemblages  with  preferred  horizontal  bedding  and  increases 
as  the  particles  become  mwe  angular.  When  the  bedding  plane  is  inclined  toward  45°  +  <j>/2,  however,  the  strength 
and  dilatioi  are  significantly  less,  particularly  when  the  particles  become  more  angular.  Assemblages  with  more 
isotropic  (or  random)  bedding  ediibit  a  mod^  increase  in  strength  and  dilation  as  particle  aspect  ratio  increases. 
Like  soU  and  other  real  granular  materials,  the  strength  and  dilation  are  lower  for  looser  numerical  assemblages. 

This  Section  attenq)ts  to  explain  wl^r  particle  sh^  and  fabric  anisotropy  influence  the  macroscopic  behavior 
using  Rowe's  stress-dilatancy  theory  for  cohesionless  granular  materials.  This  discussion  of  the  ellipse-based 
DEM  results  suggest  that  angular  particles  interlock  more  and,  as  a  result,  remold  and  dilate  significantly  less 
when  the  bedding  plane  is  inclined  with  respect  to  the  majcn*  principal  plane.  Strength  and  dilatant  behavior  is 
minimum  for  assemblages  consisting  of  angular  particles  when  the  preferred  bedding  plane  is  coaxial  with  the 
slip  planes  (45°  ±  ^/2).  As  the  bedding  plane  deviates  fixxn  the  slip  plane,  the  numerical  assemblages  dilate  more 
ai^  exhibit  more  strengdi  as  particle  intoloddng  is  mcae  easily  deteriorated  by  the  vertical  shear  stress.  Strength 
and  dilatant  behavior  is  highest  when  the  bedding  plane  is  coaxial  with  the  major  principal  plane. 


9.1  Stress-dilatancy 


Soil  and  odier  disomtinuous  materials  are  typically  anisotropic  in  nature  and  are  known  to  contract  and  dilate 
during  the  course  of  deformation.  Rowe  (1962)  used  basic  energy  principles  with  regular  and  irregular  packing 
of  two-  and  three-dimoisional  round  particles  to  develop  the  following  relationship  between  the  strength  and 
deformation  bdiavior  for  cohesionless  granular  materials  subjected  to  triaxial  compression; 


tan^(45  +^4y) 


1  +  dV/Ve^ 


(59) 


where  <|)f  is  the  internal  faction  angle  corrected  for  the  energy  due  to  expansion,  p  '  and  p  '  are  the  effective 
principal  stresses  and  dVA^e,  is  the  volumetric  strain  rate.  This  relationship  states  simply  that  the  work  done 
on  a  system  by  the  major  principal  stress  a, '  is  proportional  to  the  work  done  by  the  system  on  the  minor 
principal  stress  03'. 


Rowe's  work  shows  that  the  strength  of  granular  materials  is  related  to  the  amount  of  energy  that  is  spent 
on  dilation,  remolding  and  fiiction  as  shown  in  Figure  159.  In  order  for  an  assemblage  of  particles  to  deform, 
a  minimum  amount  of  woik  is  q)ait  oi  fiicticnal  heat.  An  additional  amount  of  energy  is  then  required  to  remold 
and  dilate  the  assemblage.  Fot  loose  assemblages,  the  additional  energy  is  spent  on  remolding  (or  particle 
rearranging)  and  on  dilating  for  dense  assemblages.  For  intermediate  assemblages,  the  additional  energy  is  spent 
on  a  crnnbinatiom  of  remolding  and  dilating. 


The  results  of  the  ellipse-based  DEM  simulations  show  that  assemblages  with  angular  particles  generally 
exhibit  hi^ier  strength  tiian  those  with  round  particles.  The  increase  in  strength  due  to  increased  particle  aspect 
ratio,  however,  is  significantly  reduced  as  the  preferred  bedding  plane  (or  interlocking  plane)  inclines  toward  the 
slip  plane  (45°  -  (|>/2  fi'om  the  major  principal  stress  q  ).  When  the  bedding  plane  and  the  slip  plane  become 
coaxial,  an  "easy"  f^ure  plane  is  formed  and  the  strength  exhibited  tends  toward  the  minimum,  particularly  for 
assemblages  consisting  of  flatter  particles.  The  DEM  results  show  that  strength  is  maximum  when  the  bedding 
plane  is  aligned  parallel  to  the  minor  principal  stress  03  and  increases  with  aspect  ratio. 
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Figure  159.  Rowe’s  (1962)  stress-dilatancy  theory  for  cohesionless  granular  materials. 

As  particle  aspect  ratio  increases,  particles  tend  to  interlock  more.  This  suggests  that  the  (|>f  curve  shown  in 
Figure  159  rises  as  particle  aspect  ratio  increases  as  more  energy  is  spent  on  remolding.  As  deformation 
continues,  assemblages  with  preferred  horizontal  bedding  (0  =  0°)  dilate  sigmficantly  more  and  achieve  higher 
values  for  <}>„«  those  with  inclined  bedding  planes. 

To  detennine  the  effect  of  increased  particle  aspect  ratio  on  the  energy  spent  on  remolding,  (|)f  was  computed 
for  the  niimwical  assemblages  using  equation  (59)  with  values  for  3**d  AVA^ej  taken  at  the  peak  shear 
stress  Qn*;.  It  was  necessary  to  use  the  values  at  because  and  (AVA^ e,)^  did  not  always  occur  at 

the  exact  moment  of  shear.  However,  Cj/aj  and  AVA^e,  were  approximately  maximum  at 

Values  obtained  fw  (t>fusing  equation  (59)  for  assemblages  with  0  and  60°  bedding  are  shown  in  Figure  160. 
The  increase  in  <|>f  due  to  particle  aspect  ratio  for  assemblages  with  0°  bedding  is  sigmficantly  higher  than  for 
those  with  60°  bedding.  The  results  indicate  that  more  energy  is  spent  on  particle  rearranging  (or  remolding) 
when  the  bedding  plane  is  aligned  parallel  to  02  (0  =  0°),  and  as  a  result,  4  3nd  4i«(  increase  sigmficantly.  A 
comparison  of  Figures  98  and  160  shows  that  both  4imx  4f  increase  by  approximately  30°  due  to  increased 
particle  aspect  ratio  fw  assanblages  with  Ixxizraital  bedding.  This  suggests  that  the  increase  in  <t)ni«x  is  due  almost 
entiiely  to  the  inoease  in  the  aiergy  required  to  rearrange  the  angular  particles  and  not  the  amount  of  energy  spent 
on  the  assemblage's  dilation  during  the  course  of  deformation. 

For  assemblages  with  60°  bedding,  however,  less  remolding  and  dilating  occur  as  the  flatter  particles  tend 
to  interlock  more  during  the  course  of  deformation  and  move  together  as  a  "rigid  bodty".  Since  less  energy  is 
spent  on  remolding  and  dilating  in  assemblages  with  60°  bedding,  the  increase  in  <j>f  4i»x  *s  much  less 
significant  than  for  those  with  preferred  horizontal  bedding. 
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Figure  160.  Effect  of  particle  aspect  ratio  (xi  <|>f  (internal  fiicticxi  angle  corrected  for  energy  spent  on  dilation) 

for  DEM  samples  with  0  and  60°  bedding. 

DEM  samples  with  raiKkxn  bedding  are  more  isotropic  and  deform  with  very  little  volume  change.  Almost 
no  volumetric  strain  is  exhibited  by  the  random  sample  consisting  of  round  particles.  However,  as  aspect  ratio 
increases,  more  interlocking  occurs  resulting  in  slightly  more  remolding  and  dilating  behavior.  The  increase  in 
the  contractant  and  dilatant  behavior  for  assemblages  with  random  bedding  is  significantly  less  than  for  those 
with  prefened  bedding.  Accordingly,  the  increase  in  <|>f  and  <!)„„„  is  also  significantly  less.  Random  samples  with 
round  particles  (a/b  =  1.0)  yield  values  for  <t»f  and  of  19.1  and  21.1°,  respectively,  while  those  with  flatter 
particles  (a/b  =  3.0)  yield  values  of  27.4  and  31.3°.  Note  that  <|>„^  for  the  random  sample  with  round  particles 
is  less  than  the  inteiparticle  sliding  fiiction  angle  <t>j,  of  26°.  This  is  due  to  the  negative  effects  of  particle  rolling 
on  the  true  fiictional  resistance. 

According  to  Rowe  and  Figure  159  looser  assemblages  dilate  less  and  remold  more  than  denser 
assemblages  resulting  in  lower  values  for  higher  values  for  ^ .  Results  of  the  ellipse-based  DEM 

simulatiois  agree  with  Rowe's  theory.  For  munerical  assemblages  with  a/b  =  2.0  and  0°  bedding,  those  that  are 
looser  yield  lower  average  values  for  (|>™„  and  slightly  higher  average  values  for  <|>f  (e^  =  0. 175,  =  40.2°,  <|>f 

=  30.9°)  than  those  that  are  denser  (eo  =  0. 158,  (j),^  =  46.7°,  (jtf  =  27.2°). 

Thoug)i  Rowe  assumes  a  sliding-type  failure,  the  true  fiiction  angle  is  generally  regarded  as  consisting 
of  both  sliding  and  rolling  fiictional  mechanisms.  It  has  been  suggested  that  particle  rolling  is  the  dominant 
deformati(xi  medianism  fcM"  assemblages  with  high  interparticle  fiiction  (Oda  et  al  1983).  The  DEM  assemblages 
with  round  particles  (a/b  =  1 .0)  and  <|)j, = 26°  yield  values  fOT  <tinw  ranging  fi’om  1 9  to  26°.  Since  <|>max  is  lower  than 
(j)^  in  some  cases,  the  difference  is  likely  due  to  increased  particle  rolling.  DEM  assemblages  with  (|»^  =  52°  yield 
higher  values  for  than  for  those  with  =  26°.  However,  the  results  show  th^|;j^  does  not  increase 
monotcmically  with  the  increase  in  which  is  in  agreement  with  Skinner  (1969).  This  lends  fiuther  proof  that 
the  negative  effects  of  particle  rolling  on  the  fiictional  resistance  and  the  overall  strength  of  granular  materials 
increases  as  the  interparticle  sliding  fiiction  increases. 
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As  mentioned  previously  in  Section  8.2.5,  the  ellipse-based  DEM  model  allows  the  decomposition  of 
contact  deformations  into  rotation  and  translation  components  (T ing  et  al  1 993).  The  incremental  quantities  are 
summed  over  every  contact  to  assess  the  relative  importance  of  rolling  and  sliding  mechanisms  within  the 
as.semhlage  Figure  161  illustrates  the  effect  of  increased  interparticle  friction  angle  on  the  amount  of  particle 
reeling  that  occurs  in  the  numoical  assemblages  during  the  course  of  deformation.  Particle  sliding  is  reduced  by 
the  inoease  in  frictional  resistance.  Consequently,  particle  rolling  must  increase  in  order  for  the  assemblage  to 
defcMm  thoeby  becoming  the  dominant  deformation  mechanism,  which  seems  to  confirm  previous  hypotheses. 


Axial  Strain  (%) 


Figure  161.  Cumulative  particle  rolling  statistics  for  DEM  samples  with  a/b  -  1.4  and  0®  bedding  with 
variations  in  the  interparticle  fricticxi  angle  4)^. 

9.2  The  Failure  Plane 

According  to  Rowe,  failure  of  discontinuous  materials  occurs  when  the  a-plane  (or  the  preferred  bedding 
plane  0)  corresponds  to  the  average  packing  geometry  of  critical  particle  groups  at  the  instant  of  collapse. 
FurthermOTe,  Rowe  states  that  the  slip  plane  is  (xiented  at  45°  -  <|>/2  with  respect  to  the  direction  of  the  major 
principal  stress  and  is  not  related  to  the  peak  strength  of  the  assemblages.  This  apparently  contradicts  the 
cemtinuum-based  Mohr-Coulcwnb  theory  prediction  of  45°  +  4>imj/2. 

Figures  162  throu^  164  plot  particle  displacement  vectors  for  assemblages  with  a/b  =  1.4  and  preferred 
bedding  0  of  0, 30  and  90°,  respectively.  The  dilacements  of  flatter  particles  (a/b  =  3.0)  with  preferred  particle 
beddingG  of 0, 30, 60  and  90°  are  shown  in  Figures  165  through  168.  Particle  displacements  are  shown  as  lines 
that  indicate  the  initial  and  final  positions  of  the  particle  centers.  Distinct  failure  planes  are  visible,  particularly 
in  assemblages  with  inclined  particle  bedding.  Since  angular  particles  interlock  more,  they  tend  to  form 
subdexnains  that  move  together  as  a  group.  When  the  preferred  bedding  plane  is  horizontal  (or  coaxial  with  the 
majex'  principal  plane),  the  interlocking  is  rapidly  deteriorated  under  the  vertical  shear  stress  as  the  particles  are 
predisposed  to  move  laterally.  As  &e  bedding  plane  inclines  toward  the  predicted  slip  plane,  however,  the 
intoioddng  does  not  deteriorate  as  easily.  Instead,  inclined  angular  particles  tend  to  remain  interlocked  and,  as 
a  result,  mwe  rigid  body-type  deformation  occurs. 
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Figure  163.  Particle 


Figure  166.  Particle  displacements  for  a  DEM  sample  with  a/b  =  3.0  and  30°  bedding  at  10%  axial  strain. 


Figure  167.  Particle  displacements  for  a  DEM  sample  with  a/b  =  3.0  and  60°  bedding  at  10%  axial  strain. 
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Figure  168. 


Figure  169. 


Particle  displacements  for  a  DEM  sample  with  a/b  =  3.0  and  90°  bedding  at  10%  axial  strain. 


07112 
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System  plot  of  DEM  sample  with  a/b  =  3.0  and  0°  bedding  prior  to  biaxial  shear. 
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Figure  170.  System  plot  of  DEM  sample  with  a/b  =  3.0  and  0®  bedding  at  10%  axial  strain. 


Figure  169  shows  a  DEM  sample  consisting  of  flat  particles  (a/b  =  3.0)  with  preferred  horizontal  bedding 
(0  =  0°)  prior  to  shear,  while  Figure  170  shows  the  same  sample  at  approximately  10%  axial  strain.  Prior  to 
shear,  the  angular  particles  interlock  and  form  an  extremely  dense  assemblage.  When  the  numerical  sample  is 
sheared,  however,  particle  interlocking  deteriorates  rapidly  under  the  applied  shear  stress  cr,  as  the  particles  are 
predisposed  to  move  horizontally  in  the  direction  of  the  minor  principal  stress  02.  Particulate  columns  form  to 
resist  O]  resulting  in  spalling-type  behavior  and  producing  more  lateral  expansion  (or  dilation).  Failure  occurs 
when  the  columns  be(X)me  unstable  and  collapse.  Similar  spalling-type  behavior  was  observed  by  Oda  et  al 
(1983)  in  biaxial  shear  tests  using  photoelastic  rubber  rods. 

At  large  axial  strains,  slip  planes  are  clearly  visible  and  possess  positive  (inclined  toward  the  right)  or 
negative  (inclined  toward  the  left)  inclination.  Particles  on  the  slip  planes  roll  and  slide  producing  large  shear 
strains,  while  those  between  the  slip  planes  remain  interlocked,  moving  together  as  a  single  particulate  Ixxty  with 
significant^  less  shear  straia  Figure  171  plots  the  contours  of  maximum  shear  strain  for  a  larger  numerical 

assemblage  consisting  of 2004  "flat"  particles  with  a/b  -  3.0  and  0°  bedding  at  approximately  18%  axial  strain. 

Figure  172  shows  a  DEM  sample  cwisisting  of  flat  particles  (a^  =  3.0)  with  preferred  60°  bedding  prior 
to  shear,  while  Figure  173  shows  the  same  san:q>le  at  ^proximatety  10%  axial  strain.  The  spalling-type  behaviOT 
observed  in  the  DEM  sample  with  0°  bedding  is  not  evident  in  the  sample  with  0  =  60°.  Assemblages  with 
inclined  bedding  dilate  less  because  itiOTe  rigid  Ixxty-type  defcamation  occurs  as  the  inclined  particles  remain 
interlocked  throughout  the  assemblage's  deformation.  Ihe  result  is  the  formation  of  a  rather  wide  single  shear 
band  in  which  the  particles  are  significantly  disturbed  due  to  rolling  and  sliding.  On  both  sides  of  the  shear  band, 
however,  very  little  remolding  aiKl  dilating  occurs  in  the  granular  wedges  as  the  particles  remain  interlocked. 


Figure  171. 


Figure  172. 


Maximum  shear  strain  contours  fcx"  DEM  sample  with  a/b  =  3.0  and  0  —  0°  at  18%  axial  strain. 


System  plot  of  DEM  sample  with  a/b  =  3.0  and  60°  bedding  prior  to  biaxial  shear. 
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Figure  173. 


Figure  174. 
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Whik  podcets  of  large  local  shear  strains  occur  in  the  sample  with  preferred  horizontal  bedding  (0  =  0®),  a 
band  of  large  shear  strains  develc^  in  the  sann)le  with  60°  bedding  as  shown  in  Figure  174.  The  maximum  shear 
strains  within  the  shear  band,  however,  are  significantly  lower  in  magnitude  (25%)  than  the  strains  experienced 
Ity  the  sample  with  0  =  0°.  It  should  be  noted,  however,  that  the  numerical  ass^blages  are  shown  at  different 
stages  of  shear  (or  axial  strains).  The  shear  strains  decrease  sigmficantly  in  the  opposing  wedges  along  the  shear 
band. 


For  isotropic  granular  materials,  the  inclination  of  the  failure  plane  m^  be  positive  or  negative  as  both 
are  equally  likely  to  occur.  The  fabric  of  the  mimical  samples  with  preferred  bedding,  however,  were  extremely 
anisotropic.  Though  positive  and  negative  inclinations  id  occur  for  the  DEM  samples  with  preferred  particle 
bedding  of  0, 30  and  90°,  the  failure  plane  inclination  was  always  negative  (inclined  toward  the  left)  fcM-  those 
with  0  =  60°. 


Particles  in  the  shear  band  roll  significantly  frOTi  their  initial  orientation,  while  those  outside  the  shear 
band  do  not  rotate  significantly  because  of  more  interlocking  and  "rigid  body"  motion.  For  a  DEM  sample  with 
a/b  =  3.0  and  60°  bedding,  particles  in  the  shear  band  that  rotate  (A0)  more  than  10°  during  the  course  of 
deformation  are  shown  in  Figure  175,  while  those  with  A0  of  more  than  20°  are  shown  in  Figure  176. 
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Figure  175.  Particles  in  a  DEM  sample  with  a/b  =  3.0  and  0  =  60°  having  a  change  in  rotation  (A0)  of 
greater  than  10°  at  10%  axial  strain 

A  nnmfiriral  assemblage  was  biaxially  sheared  in  separate  simulations  with  lower  interparticle  sliding 
fiicticxi  angles  than  26°  to  determine  whether  a  sliding-type  failure  could  be  initiated.  The  simulations  were 
p^cxmed  tnvler  the  same  conditions  except  that  <j>^  was  varied  using  values  of  5, 15  and  23°.  Fcm-  assemblages 
with  (|)^  =  5  and  15°,  the  inclinaticxi  of  the  failure  plane  was  positive  (with  the  grain)  indicating  sliding  failure, 
vriiile  the  sample  with  <})^  =  23°  initially  developed  a  positive  inclination  but  ultimately  failed  against  the  grain. 
This  seems  to  confirm  that  particle  rolling  is  the  dominant  failure  mechanism  for  assemblages  with  high  sliding 
fiiction,  even  cadsing  samples  with  preferred  60°  bedding  to  fail  against  the  grain  of  the  bedding  plane. 
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Figure  176.  Particles  in  a 


9  J  Contact  Normal  Distribution 

Results  of  the  ellipse-based  DEM  simulaticms  show  that  higher  strength  is  exhibited  by  assemblages  that 
have  initially  anisotn^ic  contact  ncxmal  distributicxis  parallel  to  the  direction  of  the  vertically  applied  shear  stress 
o,  (or  when  the  bedding  plane  is  horizontal).  Less  strength  is  exhibited  when  the  initial  distribution  anisotropy 
and  <j,  are  skewed  Fix’  assemblages  with  nicxe  isotn^icalty  distributed  contact  normals,  the  strength  is  also  less. 
During  the  course  of  deformation,  the  contact  normals  tend  to  rotate  toward  the  major  principal  stress  axis. 
Similar  fabric  evduticxi  has  been  observed  in  numerous  (nevious  physical  and  numerical  micromechanical  studies 
(Konishi  et  al  1983,  Oda  et  al  1980,  Rothenburg  and  Bathurst  1991,  Chen  and  Ishibashi  1990). 

Angular  particles  interlock  more  and  possess  more  contacts  as  indicated  by  the  increase  in  the  coordination 
numbo*.  For  horizontal  bedding  (6  =  0®),  an  increase  in  the  number  of  vertical  contact  normals  produces  a  more 
anionhryir.  fabric,  thereby  providing  more  stability  and  more  resistance  to  the  vertical  shear  stress.  For  inclined 
bedding,  however,  the  increase  in  the  number  of  contact  normals  provides  sigmficantly  less  benefit  because  the 
alignmfint  of  the  more  anisotropic  fabric  is  skewed  with  o,.  For  random  bedding,  an  increase  in  the  number  of 
particle  contacts  provides  only  slightly  more  stability  and  shear  resistance  because  the  contact  normals  are 
distributed  more  randomly  and  ewdy  a  few  more  are  aligned  toward  o,. 

Rose  diagrams  illustrating  the  contact  normal  distribution  for  a  DEM  sample  consisting  of  round  particles 
(a/b  =  1.0)  prior  to  shear  and  at  peak  strength  are  shown  in  Figure  177.  Initially,  the  contact  normals  are 
isotrc^ically  distributed  as  die  round  particles  do  not  possess  preferred  orientations.  At  peak  strength,  however, 
the  distribution  is  more  anisotrc^ic  and  (xiented  toward  the  vertical  (parallel  to  Oi).  As  the  assemblage  defexms, 
hori^/^ntal  oMitacts  are  brdeen  while  more  vertical  contacts  are  made. 


Figure  177.  Distribution  of  ccxitact  tKxmals  fix  a  DEM  sample  with  a/b  -1.0  (no  orientation)  prior  to  shear 

(left)  and  near  peak  strength  (right). 
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Figures  178  through  180  show  the  initial  and  peak  strength  contact  normal  distributions  fcx*  DEM  samples 
with  preferred  horizontal  bedding  (0  =  0®)  and  a/b  =  1.1, 1.7  and  3.0,  respectively.  Though  the  coordinaticm 
numbCT  for  the  assonblages  consisting  of  flatto’  particles  decreases  during  the  course  of  deformation,  the  majority 
of  the  lost  contacts  are  hofizcmtal.  The  maximum  fiequency  of  the  cmitact  normals  indicate  that  vertical  contact 
normals  are  still  made  even  while  the  coordination  number  decreases.  Higher  strength  is  exhibited  by 
assonbtages  consisting  of  flatto*  particles  as  the  increase  in  the  munber  of  vertical  ccmtact  normals  provide  mwe 
stability  and  shear  resistance. 


Figure  178.  Distributiai  of  contact  normals  fcx  a  DEM  sample  with  a/b  =  1. 1  and  0°  bedding  prior  to  shear 

(left)  and  near  peak  strength  (right). 


Figure  1 79.  Distribution  of  craitact  normals  fix  a  DEM  sample  with  a/b  =  1 .7  and  0°  bedding  prior  to  shear 

(left)  and  near  peak  strength  (right). 
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Figure  180.  Distributicxi  of  cxxitact  iKxmals  for  a  DEM  sample  with  a/b  -  3.0  and  0®  bedding  prior  to  shear 

(left)  and  near  peak  strength  (right). 

For  assemblages  with  preferred  60®  bedding,  the  contact  normals  are  initially  oriented  60®  from  the 
vertical  and,  as  a  result,  do  not  provide  as  much  resistance  to  a  vertical  shear  stress.  Figures  181  through  183 
showtheorientationandfrequencydistributicxi  of  contact  normals  for  DEM  samples  with  a/b  =  1.1, 1.7  and  3.0, 
respectively.  The  rose  diagrams  become  more  anisotropic  as  particle  angularity  increases.  As  the  assemblages 
are  sheared,  the  skewed  rosettes  tend  to  rotate  toward  the  vertical.  Though  the  assemblages  consisting  of  "flat" 
particles  (a/b  =  3.0)  are  more  anisotropic  and  possess  a  higher  number  of  contact  normals  than  those  consisting 
of  rounded  particles,  no  significant  increase  in  the  strength  is  exhibited  as  the  amsotropic  distribution  is  not 
parallel  to  Cj  (or  vertical). 
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Figure  181 .  DistributicKi  of  ccxitact  ixxmals  foe  a  DEM  sample  with  a^  =  1 . 1  and  60®  bedding  prior  to  shear 

(left)  and  near  peak  strength  (right). 
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Figure  182.  DistributicMi  of  contact  mxmals  for  a  DEM  sanq>le  with  a/b  =  1 .7  and  60°  bedding  prior  to  shear 

(left)  and  near  peak  strength  (right). 


Figure  183.  Distribution  of  coitact  ncxmals  fo"  a  DEM  sample  with  a/b  =  3.0  and  60°  bedding  prior  to  shear 

(left)  and  at  peak  strength  (right). 

Figures  184  and  185  show  the  initial  and  peak  strength  ccmtact  normal  distributions  for  DEM  samples 
with  30°  bedding  and  particle  aspect  ratios  a/b  of  \  A  and  3.0,  respectively.  Rose  diagrams  for  the  numerical 
assemblages  with  90°  bedding  and  a^  =  1.4  and  3.0  are  shown  in  Figures  186  and  187.  Assemblages  with  30° 
bedding  exhibit  less  strength  than  those  with  0°  bedding  due  in  part  to  the  skewed  alignment  of  the  contact 
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normals  wifli  lespect  to  a,.  The  strength  is  even  less  for  samples  with  90°  bedding  as  the  initial  horizontal  fabric 
anisotropy  jwovides  fewer  vertical  cfMtact  ncHinals  and  less  resistance  to  the  vertically  applied  o, .  As  the  samples 
are  sheared,  the  horizontal  anisotropy  deteriorates  and  re-orientates  toward  the  vertical  as  more  horizontal 
contacts  are  lost  vdiile  vertical  contacts  are  made.  For  samples  crmsisting  of  more  angular  particles  (a/b  =  3.0), 
deterioration  and  re-orientation  is  inhibited  by  the  more  anisotropic  fabric. 


Figure  184.  Distributicxi  of  contact  nomals  fix’  a  DEM  sample  with  a/b  =  1 .4  and  30°  bedding  prior  to  shear 
(left)  and  near  peak  strength  (right). 


Figure  185.  Distributicxi  of  ccxitact  normals  fix'  a  DEM  sanple  with  a/b  -  3.0  and  30°  bedding  prior  to  shear 

(left)  and  near  peak  strength  (right). 
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Figure  186.  Distribudcxi  of  coitact  nmnals  fw  a  DEM  san^le  with  a/b  =  1.4  and  90°  bedding  prior  to  shear 

(left)  and  near  peak  strength  (rig^t). 


Figure  187.  Distributicm  of  contact  normals  fw  a  DEM  san^le  with  a/b  =  3.0  and  90°  bedding  prior  to  shear 

(left)  and  near  peak  strength  (right). 

Figures  188  through  191  show  the  initial  (prior  to  shear)  and  peak  strength  contact  normal  distributions  for 
samples  with  randomly  oriented  particles  and  a/b  =  1.0,  1.1,  1.7  and  3.0,  respectively.  Though  the  random 
samples  are  much  looser  than  the  assemblages  with  preferred  bedding,  their  initial  coordination  numbers  were 
shown  to  be  very  similar.  However,  because  the  initial  fabrics  for  the  random  samples  are  significantly  less 
anisotropic,  their  evdutiai  toward  the  m^or  principal  stress  axis  occurs  more  easily.  The  fabric  evolution  occurs 
as  many  of  the  horizontal  crmtacts  are  broken  while  vertical  contacts  are  made.  Prior  to  shear  and  at  peak 
strength,  assemblages  consisting  of  flatter  particles  possess  more  vertical  and  less  horizontal  contacts  than  those 
with  rounded  particles. 
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Figure  188.  Distribution  of  OMitact  ncxmals  fix  a  DEM  sample  with  a/b  =  1 .0  (no  orientation)  prior  to  shear 
(left)  and  at  peak  strength  (right). 
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Figure  190.  Distribution  ofcontactiKxmals  for  a  DEM  sample  with  a/b=  1.7  and  rand(»n  bedding  prior  to 
shear  (left)  and  at  peak  strength  (right). 


Figure  191.  Distribudoi  of  coitact  ncxmals  for  a  DEM  sample  with  a/b  =  3.0  and  random  bedding  prior  to 

shear  (left)  and  at  peak  strength  (right). 


The  evolutirxi  of  the  contact  nomals  during  the  course  of  deformation  for  DEM  samples  with  a/b  =  2.0  and 
preferred  0  and  60°  bedding  is  shown  in  Figures  192  and  193,  respectively.  For  assemblages  with  0°  bedding, 
the  orientatio)  of  the  contact  normals  distributicm  remains  parallel  with  the  vertical  shear  stress.  For  assemblages 
with  60°  bedding,  however,  the  initial  coitact  ncxmals  are  aligned  60°  from  o,  and  as  the  sample  is  sheared,  their 
orientation  tends  to  rotate  toward  the  vertical.  Though  particle  angularity  is  the  same  for  both  assemblages,  the 
contact  normal  distribution  for  the  sample  with  0°  bedding  is  significantly  more  anisotropic  than  for  the 
assemblage  with  60°  bedding. 
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Figure  192.  Evolution  of  contact  normal  distributions  for  a  DEM  sample  with  a/b  =  2.0  and  0®  bedding 
(peak  strength  at  e,  =  4.7%). 
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Figure  193.  Evolution  of  (XHitact  normal  distributicms  for  a  DEM  sample  with  a/b  =  2.0  and  60°  bedding 
(peak  strength  at  ei  =  7.4%). 
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10.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  Report  has  presented  results  of  a  multi-year,  Air  Force  Office  of  Scientific  Research-funded  study  on 
the  mechanical  behavior  of  anisotropic  granular  systems.  A  two-dimensicmal  discrete  element  method  (DEM) 
model  using  ellipse-shaped  particles  was  used  to  investigate  the  influence  of  particle  shape  on  the  overall 
mechanical  behavior  of  granular  materials.  Numerous  particulate  assemblages  with  variations  in  particle 
angularity  and  bedding  were  subjected  to  biaxial  shear  in  numerical  simulations  under  similar  initial  conditions. 
Sc^  initial  ctHiditions  such  as  the  initial  void  ratio,  method  of  sample  formation  and  the  confining  stress  were 
varied  in  other  test  simulations  to  examine  their  effect  on  the  stress-dilatant  behavior  of  granular  assemblages. 
The  numerical  samples  oxisisted  of  dry,  weightless  and  cohesionless  ellipse-shaped  particles  that  were  generated 
in  "space".  Because  the  effects  of  gravity  arc  neglected,  the  significance  of  the  orientation  of  the  preferred 
bedding  plane  is  with  respect  to  the  laterally-applied  principal  (or  confining)  stresses. 

The  comparisons  of  the  numerical  test  results  to  those  from  physical  tests  on  sand  and  other  granular 
show  that  the  ellipse-based  DEM  model  can  realistically  simulate  real  soil  behavior.  Good  agreement 
was  found  with  the  circular  rod  data  of  Chapuis  (1976).  Thou|^  good  quantitative  agreement  between  the  data 
reported  by  Konishi  et  al  (1983)  and  that  yielded  by  the  ellipse-based  DEM  model  was  not  achieved,  the 
differences  are  aittributed  to  the  subtle  but  significant  differences  in  the  shape  of  the  particles  used  in  both 
micnxnedianical  investigaticms.  The  ellipse-based  DEM  model  used  in  the  current  numerical  study  yields  stress- 
dilatant  behavior  results  that  are  quantitatively  closer  to  real  soil  and  other  granular  materials  than  those  yielded 
by  other  DEM  models  using  round  or  elliptic  particle  shapes. 

The  results  of  the  biaxial  shear  simulations  using  the  ellipse-based  DEM  model  show  that  the  strength  and 
dilatant  bdiavior  of  granular  materials  generally  increases  as  the  particulate  constituents  become  more  angular. 
Upo)  deposititxi,  angular  particles  mty  beoxne  aligned  producing  dense  assemblages  with  extremely  anisotropic 
fabrics  or  they  may  be  randcMtily  distributed  producing  loose  assemblages  with  more  isotropic  fabrics.  Results 
yielded  by  the  ellipse-based  DEM  model's  biaxial  shear  simulations  show  that  fabric  anisotropy  significantly 
influences  (he  strength  and  (Mcxmation  behavior  exhibited  by  the  particulate  assemblages.  Since  round  particles 
do  not  possess  (mentation,  the  fabric  tends  to  be  mcx’e  isotropic.  Fabric  anisotropy  increases  as  particle  aspect 
ratio  increases. 

The  numerical  simulati(m  results  show  that  the  strength  of  granular  materials  is  lowest  for  assemblages 
consisting  of  round  particles  and  increases  as  particle  angularity  increases.  For  assemblages  consisting  of 
particles  with  similar  shapes,  the  strength  and  dlatant  behavior  exhibited  is  highest  when  particle  bedding  is 
normal  to  the  iqiplied  shear  stress  (or  horizontal  for  the  current  study).  The  exhibited  strength  and  dilatant 
behavior  decreases  as  the  bedding  plane  inclines  toward  the  predicted  failure  plane  defined  by  Rowe's  stress- 
dilatancy  thecny  (45®  -  <|>/2  from  applied  shear  stress  cn ).  For  assemblages  consisting  of  flatter  particles,  the 
strength  and  dilatant  behavicn^  are  lowest  when  the  bedding  plane  and  45®  -  <^/2  are  coaxial.  The  results  of  the 
numerical  simulations  show  that  stress-dilatancy  behavior  of  granular  materials  increases  as  the  bedding  (or 
interkxking)  plane  deviates  from  the  predicted  slip  plane. 

As  particle  angularity  increases,  mote  interioddng  occurs.  Ckxisequently,  the  numerical  granular  assemblages 
experience  more  remolding  and  longer  peri(xls  of  initial  contraction.  After  the  extended  initial  contractant 
behavior,  assonblages  with  preferred  horizontal  bedding  (0  =  0®)  ultimately  dilate  more  achieving  higher 
volumetric  strains  than  assemblages  consisting  of  more  rounded  particles.  As  a  result,  assemblages  with 
preferred  horizcmtal  bedding  exhibit  higher  strength  as  more  energy  is  required  to  deform  the  assemblage  by 
expansion. 

Though  the  assemblages  with  preferred  60°  bedding  are  extremely  dense,  they  exhibit  behavior  typically 
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^ggnriatPiH  with  loosc  assemblages.  Because  the  bedding  plane  of  assemblages  with  preferred  60“  bedding  arc 
^prmdmatety  coaxial  with  the  predicted  slip  plane  (45“  -  failure  is  facilitated  more  easily.  Consequently, 
the  assemblages'  deformation  more  closely  resembles  rigid  body  motion  of  opposing  particulate  wedges  along 
the  slip  plane  that  develops  as  a  distinct  shear  band.  Strength  and  dilatant  behaviors  for  these  assemblages 
decrease  with  innipAgfri  particle  angularity  as  less  aiergy  is  spait  cm  expansion  during  the  course  of  shear-induced 
defcxmatioa  Results  of  the  numoical  simulations  show  that  Icxiser  assemblages  having  isotropic  or  anisotropic 
orientation  fabrics  exhibit  less  strength  and  dilatant  behavior  than  denser  assemblages  as  more  energy  is  spent 
cm  remolding  than  cm  dilatiem. 

Assonblages  with  piefened  hcxizimtal  bedding  (0  =  0“)  that  arc  cemfined  laterally  by  higher  stresses  require 
more  energy  to  expand  against  the  increased  confining  stress  and  as  a  result  exhibit  higher  strength.  The 
in<^rftngflH  rrmfining  stTcss,  howevcT,  provides  more  resistance  to  lateral  expansion  resulting  in  less  dilatant 
behavior. 

For  granular  assemblages  with  low  interparticle  friction,  results  of  the  numerical  simulations  show  that 
particle  sliHing  is  the  dominant  deformation  mechanism.  However,  the  numerical  simulations  show  that  particle 
rolling  is  the  dominant  deformation  mechanism  in  assemblages  with  high  interparticle  frictiem.  Consequently, 
as  the  interparticle  friction  increases,  the  strength  does  not  increase  monotonically.  As  the  interparticle  friction 
increases,  assemblages  yield  peak  macroscopic  internal  friction  angles  that  are  lower  than  their  interparticle 
sliding  friction  angles. 

angular  particles  possess  more  interlocking,  the  negative  effects  of  particle  rolling  on  the  strength 
of  the  assemblage  are  significantly  reduced.  For  granular  assemblages  consisting  of  round  particles,  particle 
motion  due  to  rotation  is  about  2.5  times  larger  than  particle  translation.  As  particle  angularity  increases,  the 
relative  importance  of  particle  motion  due  to  translation  and  rotation  become  virtually  equal.  Rotation  and 
translation  statistics  are  similar  for  assemblages  with  both  preferred  0  (horizontal)  and  30°  bedding  suggesting 
that  the  relative  importance  of  rotation  and  translation  for  assemblages  with  variations  in  particle  angularity  is 
independent  of  the  orientation  of  the  preferred  bedding  plane. 

Results  of  the  biaxial  shear  simulatkxis  yielded  Ity  the  ellipse-based  DEM  model  show  that  assemblages  with 
preferred  bedding  possess  more  anisotropic  distribution  of  contact  normals.  As  particle  angularity  increases, 
interlocking  increases  producing  a  more  anisotropic  fabric.  During  the  course  of  deformation,  the  distribution 
of  contact  normal  orientations  tend  to  rotate  toward  the  vertical  (or  the  direction  of  the  applied  shear  stress) 
providing  incfeased  shear  resistance  as  horizontal  contacts  are  broken  while  vertical  contacts  are  made. 

Anisotropic  fabrics  due  to  increased  particle  interlocking  are  deteriorated  more  rapidly  in  assemblages  with 
preferred  0°  bedding  as  the  particles  are  predisposed  to  move  laterally  since  they  possess  horizontal  major  axes 
that  are  coaxial  with  the  horizontally-applied  confining  stress.  Consequently,  more  dilatant  and  spalling-type 
behavior  occurs  resulting  in  an  increase  in  the  ability  to  resist  a  vertically-applied  shear  stress.  Because  of  the 
mrrwflgfd  dilatant  behavior  of  assemblages  with  preferred  horizontal  particle  bedding,  the  coordination  number 
decreases  significantly  during  the  course  of  defamation. 

Shear  resistance  is  significantly  less  as  the  {xeferred  bedding  plane  inclines  toward  45°  -  <t>/2  (measured  from 
<y,).  Particulate  columns  that  develop  to  provide  shear  resistance  are  considerably  less  stable  than  those  with 
preferred  0°  bedding.  Particle  interlocking  increases  during  the  course  of  deformation  resulting  in  less  energy 
required  to  laterally  expand  assemblages  with  iiKlined  bedding.  Cwisequently,  the  coordination  number  increases 
as  the  granular  assemblages  deform. 
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The  ellipse-based  DEM  model  used  for  the  current  stucty  has  been  shown  to  be  an  effective  tool  for 
realistically  modeling  the  microstructure  and  overall  mechanical  behavior  of  granular  materials.  Because 
discrepancies  exist  in  the  comparison  of  the  results  reported  by  Konishi  et  al  (1983)  and  those  yielded  by  the 
ellipse-based  DEM  model,  it  is  recommended  that  the  model's  validation  continue  to  include  data  comparisons 
widi  physical  tests  that  specifically  use  ellipse-shaped  particles.  Such  a  study  is  currently  underwt^  using  the 
results  of  interfacial  shear  tests  between  a  rough  inteiface  an  assemblage  of  carefully  machined  circular  and 
elliptical  rods. 

Gradation  of  the  numerical  assemblages  used  in  the  current  investigation  was  limited  to  an  equal  number  of 
three  different-sized  particles  having  the  same  shape  (or  aspect  ratio).  Though  variations  in  the  gradation  of 
particulate  materials  can  be  limitless,  it  is  recommended  that  more  variation  in  the  assemblage's  gradation  be 
investigated.  Variations  in  gradation  should  include  mixing  of  particle  shapes  in  addition  to  finer  and  coarser 
distributicMis  of  particle  sizes. 

Additional  woik  is  progtessing  to  investigate  the  effect  of  boundaries  (rigid,  flexible,  periodic).  As  well  the 
ellipse-based  DEM  is  being  applied  to  a  limited  number  of  soil-structure  interaction  problems  in  geotechnical 
engineering.  Other  efforts  are  currently  being  made  to  enable  the  model  to  simulate  particle-fluid  interactions 
and  multiphase  materials  such  as  fiozen  soil  systems. 

The  development  of  the  ellipse-based  DEM  was  funded  initially  by  the  Natural  Sciences  and  Engineering 
Research  Council  of  Canada  through  an  operating  grant  to  the  Principal  Investigator  while  at  the  University  of 
Taooto.  Sipport  for  the  woric  described  in  the  Report  was  provided  by  the  U.S.  Air  Force  Office  of  Scientific 
Research,  Aerospace  ScioKes  Directorate  through  grant  F49620-92-J0017.  All  sources  of  funding  and  support 
are  very  gratefully  acknowledged. 
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APPENDIX  II.  SYMBOLS 
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a„  Uj,  aj 
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''v 
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major  and  minor  axes  for  ellipse 

normal  and  shear  contact  viscous  damping  coefficients 

coefficients  of  general  ellipse  function 

center  of  ellipse  in  wc»’ld  coordinates 

coefficient  of  restitution 

unit  normal,  outward  fm*  ellipses,  inward  for  wall 
component  of  force  fcx*  contact  c 
nmmal  and  tangential  forces  at  a  contact 
polar  momrat  of  inertia  of  particle 
normal  and  tangoitial  contact  spring  stiffiiesses 

/*  component  of  contact  vectw  (particle  centroid  -  contact  location)  for  contact  c 
mass  of  particle 

slope  and  j^-intercept  of  wall  in  homogeneous  local  coordinates  of  ellipse  / 

moment  about  point  o 

total  number  of  particles  in  a  region 

constants  in  ellipse-ellipse  intersection  equation 

sum  of  squares  of  errcn? 

unit  tangential  vector  to  ellipse 

maximum  contact  normal  and  shear  strengths 

displacement  components  of  /'*  particle 

volume 

unit  vector  along  wall  from  endpoint  1  to  endpoint  2 

Cartesian  coordinate  variables 

endpoints  of  wall 

Coordinates  of  contact 

homogeneous  coordinate  variables 

coefficients  of  displacement  field  m  =  «)  +  _y, 

coefficients  of  displacement  field  v=J3,+fi2  +  fij  yi 

average  strain  tensor 

average  stress  tensor 

rotation  with  respect  to  world  coordinate  systmn 


159 


